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ABSTRACT 

We substantially update the capabilities of the open source software package Modules for Experiments in 
Stellar Astrophysics (MESA), and its one-dimensional stellar evolution module, MESA star. Improvements in 
MESA star's ability to model the evolution of giant planets now extends its applicability down to masses as low 
as one-tenth that of Jupiter The dramatic improvement in asteroseismology enabled by the space-based Kepler 
and CoRoT missions motivates our full coupling of the ADIPLS adiabatic pulsation code with MESA star. This 
also motivates a numerical recasting of the Ledoux criterion that is more easily implemented when many nuclei 
are present at non-negligible abundances. This impacts the way in which MESA star calculates semi-convective 
and thermohaline mixing. We exhibit the evolution of 3-8 M© stars through the end of core He burning, 
the onset of He thermal pulses, and arrival on the white dwarf cooling sequence. We implement diffusion 
of angular momentum and chemical abundances that enable calculations of rotating-star models, which we 
compare thoroughly with earlier work. We introduce a new treatment of radiation-dominated envelopes that 
allows the uninterrupted evolution of massive stars to core collapse. This enables the generation of new sets of 
supernovae, long gamma-ray burst, and pair-instability progenitor models. We substantially modify the way in 
which MESA star solves the fully coupled stellar structure and composition equations, and we show how this 
has improved MESA's performance scaling on multi-core processors. Updates to the modules for equation of 
state, opacity, nuclear reaction rates, and atmospheric boundary conditions are also provided. We describe the 
MESA Software Development Kit (SDK) that packages all the required components needed to form a unified 
and maintained build environment for MESA. We also highlight a few tools developed by the community for 
rapid visualization of MESA star results. 
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evolution — stars: massive — stars: rotation 
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1. INTRODUCTION 

As the most commonly observed objects, stars remain at 
the forefront of astrophysical research. Advances in opti- 
cal detector technology, computer processing power, and data 
storage capability have enabled new sky surveys (e.g., the 
Sloan Digital Sky Survey; York et al. 2000); triggered many 
new optical transient surveys, such as the Palomar Transient 
Factory (Law et al. 2009) and Pan-STARRSl (Kaiser et al. 
2010); and allowed for space missions (e.g., Kepler; Koch 
et al. 2010) that continuously monitor more than 100,000 
stars. The stellar discoveries from these surveys include rev- 
elations about rare stars, unusual explosive outcomes, and 
remarkably complex binaries. The immediate future holds 
tremendous promise, as both the space-based survey Gaia (de 
Bruijne 2012; Liu et al. 2012) and the ground based Large 
Synoptic Survey Telescope (LSST; Ivezic et al. 2008) come 
to fruition. 

These developments have created a new demand for a re- 
liable and publicly available research and education tool in 
computational stellar astrophysics. We introduced the open 
source community tool MESA (Paxton et al. 2010, hereafter 
Paper I) to meet these new demands. This first "instrument" 
paper described the design and implementation of MESA mod- 
ules for numerics, microphysics, and macrophysics, and intro- 
duced the stellar evolution module, MESA star. We presented 
a multitude of tests and code comparisons that served as our 
initial verification and demonstrated MESA star's initial capa- 
bilities. Since Paper I, MESA has attracted over 500 registered 
users, witnessed over 5,000 downloads from http : //mesa . 
sourceforge.net/, started an annual Summer School pro- 
gram, and provided a portal (http : //mesastar . org) for the 
community to openly share knowledge (e.g., the specific set- 
tings for a published MESA star run), codes, and publications. 

This paper describes the major new MESA capabilities for 
modeling giant planets, asteroseismology, and the treatment 
of rotation and evolution of massive stars. We also describe 



numerous advances since Paper I. These include the incorpo- 
ration of composition gradients in the determination of con- 
vective mixing and additional verification for evolution of in- 
termediate mass stars and the white dwarfs they create. 

Our improvements to MESA star for gas giant planets were 
motivated by the dramatic growth in this field. Over 800 ex- 
oplanets have been confirmed, and their study has prompted 
enormous progress in our understanding of the formation and 
migration of giant planets, and of the importance of factors 
such as stellar mass (Laughlin et al. 2004; Alibert et al. 201 1; 
Boss 2011), composition (Fischer & Valenti 2005; Young 
et al. 2012), and binarity (Patience et al. 2002; Mugrauer & 
Neuhauser 2009; Roell et al. 2012). Puzzles remain, though, 
both in our solar system and in the studies of the plethora 
of these newly discovered exoplanets, including the charac- 
teristics of the planet-hosting stars and the interiors, atmo- 
spheres, surface gravities, temperatures, and compositions of 
the planets (e.g., Udry & Santos 2007; Seager & Deming 
2010). Many of these variations can now be numerically ex- 
plored, as can the incorporation of an inert core in an other- 
wise regular gas giant and the impact of irradiation. 

The ability to infer stellar properties (e.g., mass, radius, in- 
ternal state, and rotation) from measurements of the radial and 
non-radial oscillation modes has been dramatically improved 
by two space-based optical telescopes (Convection Rotation 
and Planetary Transits, CoRoT; Baglin et al. 2009 and Ke- 
pler, Borucki et al. 2009). The high cadences and precision 
(often better than ten parts per million) reveal and accurately 
measure multitudes of oscillation frequencies for over 10,000 
stars, substantially raising the need for accurate and efficient 
computations of stellar mode frequencies and the resulting 
eigenfunctions. The intrinsic flexibility of MESA star allows 
for the exploration of model-space required to precisely infer 
stellar properties from the observed frequencies. 

An important new addition to MESA is the incorporation of 
stellar rotation and magnetic fields in radiative regions. As 
stars are not solid bodies, they undergo radial differential rota- 
tion (Thompson et al. 2003; Balbus et al. 2012) and also rotate 
at different angular velocities at different latitudes (Ruediger 
et al. 1998; Bonanno et al. 2007; Kiiker et al. 2011). These 
rotational shears have a significant impact on the evolution 
of the stellar magnetic field. Despite the resulting 3D nature 
of magnetism and rotation, the stellar evolution community 
has come a long way in understanding stars with ID simu- 
lations (Meynet & Maeder 1997; Langer et al. 1999; Maeder 
& Meynet 2000b; Heger & Langer 2000; Hirschi et al. 2004; 
Cantiello & Langer 2010), thus motivating our need to fully 
incorporate rotation within MESA. The new flexibility in an- 
gular momentum transport mechanisms allows for numerical 
exploration of alternate rotational outcomes should the obser- 
vations (e.g., asteroseismology) require it. 

The paper is outlined as follows. Section 2 describes the 
new capability of MESA to evolve models of giant planets, 
while §3 discusses the new asteroseismology capabilities. 
The MESA implementation of composition gradients in stellar 
interiors and their impact on convective mixing is described in 
§4. The status of the evolution of intermediate mass stars and 
the MESA star construction and evolution of white dwarfs is 
described in §5. The new capabilities for evolving rotating 
stars is described in §6. The onset of near Eddington lumi- 
nosities and radiation pressure dominance in the envelopes of 
evolving massive stars has been a challenge for many stellar 
evolution codes ever since the realization of the iron opac- 
ity bump at logT * 5.3 (Iglesias et al. 1992). We discuss 
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in §7 the resulting improvements for evolving massive stars. 
This allows for the uninterrupted evolution of rotating massive 
stars to the onset of core collapse. We conclude in §8 by high- 
lighting where additional improvements to MESA are likely to 
occur in the near future. Appendix A describes the many im- 
provements to the physics modules since Paper I; Appendix B 
presents "nuts and bolts" information on the primary compo- 
nents of evolution calculations; and Appendix C presents the 
MESA Softwai-e Development Kit (SDK). All of our symbols 
are defined in Table 1. We denote components of MESA, such 
as modules and routines, in Courier font, e.g., evolve_star. 

Table 1 

Variable Index. 



Table 1 
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Name 


Description 


First Appears 


A 


atomic mass number 


Section A.l 


A; 


mass excess of the (tli isotope 


Section A.l 


V 


wind mass loss coefficient 


Section 5.2 




day-Side nux incident on an irradiated planet 


Section 2.4 


r 


Coulomb coupling parameter 


Section 5.2 


i 


specific angular momentum 


Section B.6 


K 


opacity 


Section 2.1 


L 


stellar luminosity 


Section 3.2 


in 


Lagrangian mass coordinate 


Section 2.1 


M 


stellar mass 


Section 2.1 


N 


Brunt- Vaisala frequency 


Section 3.3 


"i 


number density of the r'th isotope 


Section A.l 


y 


turbulent viscosity 


Section B.6 


r 


radial coordinate 


Section 2.4 


R 


total stellar radius 


Section 2.1 


P 


baryon mass density 


Section A.l 


S 


specific entropy 


Section 2.1 




mass column 


Section 2.4 




depth for heating from irradiation 


Section 2.4 


r 


optical depth 


Section 5.2 


Vc 


magnitude of changes during a timestep 


Section B.3 


Vf 


target value for 


Section B.i 


W 


atomic weight 


Section A. 1 


X 


H mass fraction 


Section 3.2 


Xi 


baryon mass fraction of the /'th isotope 


Section 3.3 


/ 


He mass fraction 


Section 2 


Ye 


electrons per baryon (Z/A) 


Section A. 1 


Yi 


abundance of the rth isotope 


Section A. 1 


z 


metallicity 


Section 2 


Z 


atomic number 


Section A. 1 


ttMLT 


mixing length parameter 


Section 2.2 


fsc 


semiconvection efficiency parameter 


Section 4.1 


ffth 


thermohaline efficiency parameter 


Section 4.2 


«V 


smoothing parameter for MLT++ 


Section 7.2 


ffV 


MLT++ parameter used in construction ot ffy 


Section 7.2 




min(P/Pg,,) 


Section 7.2 


Xp 


(din P/d In p)t 


Section 3.3 


Xt 


(d\nP/dln T)p 


Section 3.3 


Cp 


specific heat at constant pressure 


Section 4.1 


c\ 


adiabatic sound speed 


Section 3 


Av 


large frequency separation of pulsation modes 


Section 3.2 


£>ov 


overshoot dift'usion coefficient 


Section 3.3 


Ou, 


thermohaline diff'usion coefficient 


Section 4.2 


£fx 


Fermi energy at center 


Section 2.2 


%rav 


gravitational heating rate 


Section 5.2 


^nuc 


nuclear heating rate 


Section A.4 


^conv 


convective flux 


Section 7.1 


fov 


convective overshoot parameter 


Section 3.2 


frad 


radiative flux 


Section 7.1 


/v 


reduction factor for i5v 


Section 7.2 


Tl 


{din P/d In p)s 


Section 3 




opacity for thermal radiation orig. in planet 


Section 2.4 


Ky 


opacity for irradiation from star 


Section 2.4 


kB 


Boltzmann constant 


Section 2.2 




accretion luminosity 


Section 5.3 


■^max 


max(L,-ad/LEdd) 


Section 7.2 


Lc 


core luminosity 


Section 2.3 


^Edd 


Eddington Luminosity 


Section 6.3 


iogg 


log surface gravity 


Section A. 5 


^onset 


lum. at onset of convection 


Section 7.1 



Name 


Description 


First Appears 


Lrad 


radiative luminosity 


Section 7.1 


^inv 


lum. at which density inversion occurs 


Section 7.1 


mu 


atomic mass unit 


Section A. 1 


Mc 


core mass 


Section 2.3 


M 


mass-loss rate 


Section 5.3 


Mm 


modeled mass 


Section B.4 


A'a 


Avogadro number 


Section 2.2 


Vad 


adiabatic temperature gradient 


Section 3.3 


Vl 


Ledoux criterion 


Section 4.1 


Vrad 


radiative temperature gradient 


Section 3.3 


Vr 


actual temperature gradient 


Section 3.3 




baryon density 


Section A. 1 




frequency of maximum power 


Section 3.2 


n" 


surface angular velocity 


Section 6.3 


to 


angular velocity 


Section 6. 1 




surface critical angular velocity 


Section 6.3 


Pc 


central pressure 


Section 2. 1 


^as 


gas pressure 


Section 5.3 


flad 


radiation pressure 


Section 5.3 


Rc 


core radius 


Section 2.3 


Pc 


central density 


Section 2.1 


Ap 


pressure scale height 


Section 3.3 


CSB 


Stefan-Boltzmann constant 


Section 2. 1 


Sr 


Lamb frequency 


Section 3 


av 


superadiabaticity, V7- - V^d 


Section 7.2 


l5v.thrcsh 


controls when MLT++ is applied 


Section 7.2 


Tc 


central temperature 


Section 2.1 


Teff 


effective temperature 


Section 2.1 


St 


numerical timestep 


Section 5.3 


tkh 


thermal (Kelvin-Helmholtz) timescale 


Section 6.3 




equatorial velocity 


Section 6 



2. GIANT PLANETS AND LOW-MASS STARS 

Evolutionary models of giant planets and low-mass stars 
differ from their higher-mass stellar counterparts in both the 
microphysics needed to describe the interior and the role of 
stellar irradiation in the outer boundary condition. For masses 
M < 84 Mj, hydrogen burning is insufficient to prevent cool- 
ing and contraction. Deuterium burning can briefly slow the 
cooUng for M > 13 Mj, where Mj = 9.54 x 10""* Mq is 
Jupiter's mass, but has a negligible influence on the cooling 
for smaller masses. Hence nuclear burning can be ignored in 
the planetary mass regime. 

For hydrogen-helium rich objects with M » Mj, an ideal 
gas equation of state (EOS), with arbitrary degeneracy, is a 
good approximation while for M < Mj particle interactions 
play an important role. Specifically, pressure ionization of 
hydrogen at p ^ 1 g cm and r ^ 10"* K causes a sudden 
change from a H2-dominated phase to an ionized phase. MESA 
employs the Saumon et al. (1995) equation of state (SCVH 
EOS), smoothly interpolated from the low to high pressure 
phase, for this complicated region of parameter space where 
thermal, Fermi, and electrostatic energies may all be compara- 
ble. At the low temperatures in planetary atmospheres, abun- 
dant species such as CNO atoms will be in molecular form, 
and may condense into clouds. MESA does not follow the tran- 
sition from atomic to molecular form for these species in the 
EOS — they are currently included by increasing the helium 
abundance from Y to Y + Z when calling the SCVH EOS. 
MESA does, however, include the effect of molecules in the 
Rosseland opacities. Currently, the Ferguson et al. (2005) and 
Freedman et al. (2008) tables, which include the opacity from 
molecules, but ignore condensates, are available. 

Lastly, for planets in close-in orbits about their parent star, 
the external irradiation flux may be orders of magnitude larger 
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than the cooling flux from the planet's interior. This may dra- 
matically increase the surface temperature and aff'ect the outer 
boundary condition. MESA now implements several options 
for this surface heating, including the flexibility to include 
user-supplied prescriptions. 

In the following subsections, we discuss a new MESA mod- 
ule that creates initial models in the planetary mass range 
M ^ 0.1-lOMj, and present a suite of evolutionary calcu- 
lations. We discuss how surface irradiation may be included, 
as well as an inert core at the center of the planet. We also 
show what MESA star yields for the mass-radius relation for 
sub-solar mass stars in §2.5. 

2.1. Construction of Starting Models 

For stellar mass objects, the prejnsjnodel routine con- 
structs pre-main-sequence (PMS) models assuming L(r) oc m, 
where L(r) is the luminosity at radius r, by iterating on the 
starting conditions at the center to find a model with a given 
M and central temperature 7^. This PMS routine works well 
for M > 0.03 Mq, but lower masses may not converge when 
the guess for central density pc and luminosity are not close 
enough to the (unknown) true values. As a result, it is diffi- 
cult and time consuming to create models with M < 0.03 M© 
using the same routine for giant planets as for stars. 

A new routine called create_initial_model builds a 
model of given M and radius R using an adiabatic temperature 
profile. Given the central pressure and specific entropy S , 
the equation of hydrostatic balance is integrated outward, and 
the temperature at each step determined from the equation of 
state using T - T(P,S). The values of and S are iterated 
to attain the desired M and R. The luminosity profile is then 
derived treating S as constant in space for the fully convective 
planet (e.g., Ushomirsky et al. 1998), so 
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dm'T(m') = -L(m). (1) 



The luminosity at the surface, L(M), is estimated using the 
radius R and temperature Tgff at the r = KP/g -2/3 point as 
L(M) - AnR^cTsBT^ff Given L(M), the luminosity at interior 
points is found by 



L{m) = L(M) 



J" dm' T(m') 
j^^ dm' T(m')^ 



This procedure works well for M down to 
a range of initial radii. 



(2) 



0. 1 Mj and over 



2.2. Evolutionary Calculations 

Figures 1 and 2 show evolutionary calculations for 
models with masses M - 0.2-20 Mj. All models 
were evolved for 20 Gyr. The initial models from 
create_initial_model had a large radius R = 5Rj, where 
Rj - 7.192 X 10^ cm is the equatorial radius of Jupiter. 
The other parameters used are Y = 0.27, Z - 0.02 
and aMLT = 2. The opacity and EOS tables used are 
eos_£ile_pref ix - mesa, kappa_f ile_pref ix = gs98 and 
kappa_lowT_pre£ix = lowT_Freedmanll. The atmosphere 
model is which_atm_option - simple.photosphere. 

Figure 1 is a low mass extension of Figure 16 from Pa- 
per I, showing evolution in the pc-T^ plane. Each track (solid 
black curve) is labeled on the left by the planet's mass, and 
evolution goes from left to right. Initially the planet is non- 
degenerate and contraction increases both pc R^^ and 
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Figure 1. The solid blaclc lines show versus pc during the evolution. 
Each line is labeled on the left by the mass in units of Mj . The dotted red 
lines show constant values of log(age[yr]), labeled at the base of each line. 
The blue dashed lines show fixed values of 5 / (Wa^b ), labeled at the top of 
each line. The large black dots show the position of maximum along the 
evolutionary track. 




M/M,, 



Figure 2. Radius versus mass iso-contours from a suite of evolutionary cal- 
culations. The solid red lines show /?/Rj versus M/Mj at fixed values of 
log(age[yr]), labeled on the left of each curve. The dashed blue curves are 
for fixed entropy, with each curve labeled by 5/(A'a^b) on the right. The 
dotted black curves are for fixed luminosity, with each curve labeled by 
log(L[erg s"']) above M = 1 Mj. The green curve at the bottom is the 
7 = M-R relation from Zapolsky & Salpeter (1969) for a solar mixture of 
H and He. 
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Tc oc oc p^^ . A maximum % is reached when ^bTL ~ £^f,c, 
where /ip.c is the electron Fermi energy at the center, beyond 
which pc approaches a constant as decreases further Ignor- 
ing Coulomb interactions in the EOS, 5 is a function of the 
electron degeneracy parameter fif-jk^T, where jUg is the elec- 
tron chemical potential and all models should have maximum 
k^T^ ~ Ep.c at the same S . The line labeled 5/(A^a^b) = 10.3 
indeed coincides with maximum down to M ^ 1 Mj, but 
at smaller masses where non-ideal effects are more important, 
maximum occurs when 5/(A^fcB) < 10.3. Also shown in 
Figure 1 are lines of constant age, shown as dotted red lines, 
and labeled on the bottom of the plot. 

The same evolutionary calculations are used in Figure 2 to 
show radius versus mass at fixed values of age, entropy or lu- 
minosity. At late times, or low entropy and luminosity, the 
radius approaches the zero-temperature value (green curve; 
Zapolsky & Salpeter 1969) for which thermal support is in- 
significant. The maximum radius occurs where gravitational 
and Coulomb energies, per ion, are comparable. The solid 
red lines, labeled by age on the left, show that contraction 
down to R ^ 1.5 Rj is rapid, taking less than 10 Myr for 
M < lOMj. This initial rapid cooling phase occurs because 
the initial luminosity is orders of magnitude higher than the 
luminosity around one Gyr. This can been seen in the black 
dotted contours of constant log(L), where L is larger by a fac- 
tor of 100 for = 1.3 Rj and lO'* for R ^ 1.7 Rj, as compared 
to R - 1.1 Rj. The blue dashed lines show contours of con- 
stant entropy, labeled on the right by ^/(A^a^b). 

2.3. Implementation of Inert Cores 

In the core accretion model of planet formation (e.g.. Pol- 
lack et al. 1996; Hubickyj et al. 2005), a rock/ice core is first 
assembled. Once this core grows to ~ 10 Earth masses, it can 
initiate rapid accretion of nebular gas, which could then dom- 
inate the mass of the planet. For studies of planetary radii, 
a central core composed of high mean molecular weight ma- 
terial can decrease the radius of the planet by a significant 
amount (^ 0. 1-0.2 Rj). The MESA star inert core feature 
allows one to add a core of specified mass and radius 
7?c, or more conveniently, density p^. A luminosity may 
also be specified, although the high mean molecular weight 
of the core, as compared to the overlying H/He envelope, im- 
plies that even large cores will tend to have small heat content 
(Fortney et al. 2006). This inert core is not presently evolved 
in any way, and changes in during evolution are neglected 
as R changes. 

2.4. Irradiation 

Surface heating by stellar irradiation changes the boundary 
condition for the planet's cooling and contraction. This mod- 
ifies the planetary radius versus age for exoplanets at orbital 
separation < 0.1 AU. MESA provides several ways to imple- 
ment surface heating with varying degrees of fidelity to the 
true solution. These presently include: 

a) An energy generation rate e - F*/4Sv, applied 
in the outer mass column E < E*. Here is the 

rR 

day-side flux from the star, and S(r) = dr p(r ) 
is the mass column. In steady-state, this gener- 
ates an outward flux Fi,/4, which is meant to sim- 
ulate the angle-averaged flux over the planetary sur- 
face. The parameters F^, and are specified through 
the user-specified variables irradiation_flux and 



2 
1.9 
1.8 



1.7 



1.4 
1.3 
1.2 
1.1 



0.01 0.1 1 10 

age(Gyr) 

Figure 3. Radius versus age for the planet HD 209458b. The solid red lines 
are for MESA, using the grey_irradiated atmosphere model. The dotted 
black Hnes show the CEPAM code results. The dashed blue lines show the 
MESA calculation using the F^,-'Z^, surface heat source. The data point with 
error bars is the observed value of the radius for HD 209458b quoted in Guil- 
lot (2010). The two sets of curves are deep heating (upper three curves) and 
shallow heating (lower three curves). 

coluinn_depth_for_irradiation, making this the 
simplest method to use. This heating mechanism rep- 
resents absorption of stellar optical radiation well be- 
low the photosphere of the planet's thermal radiation 
and gives rise to greenhouse heating of the atmosphere 
where e 0. 

b) mesa's grey_irradiated atmosphere model (see 
also §A.5) implements the angle-averaged temperature 
profile of Guillot (2010). This approximate solution to 
the transfer equation assumes two frequency bands: op- 
tical radiation from the star (with user-specified opacity 
Ky) and thermal radiation originating in the planet (with 
user-specified opacity /Cth). The temperature profile is 
derived using the Eddington approximation, assuming 
an external flux from the star as well as a flux from the 
planetary interior. This is the only MESA atmosphere 
model that uses pressure instead of optical depth to de- 
termine the surface boundary condition. As this pres- 
sure may be relatively deep in the atmosphere, a connec- 
tion to the radius may be required to give either the ver- 
tical thermal photosphere, or the optical photosphere in 
transit along a chord. Lastly, the relax_irradiation 
routine improves initial convergence by providing a 
starting model closer to the irradiated one. 

c) Finally, MESA allows user-specified heating functions 
(e.g., Fi,-l,i, surface heating) or atmosphere models 
(e.g., grey_irradiated). User-supplied routines may 
be easily implemented by using the other_energy 
module. 

Figure 3 shows radius versus age for the planet HD 
209458b (Guillot 2010). The two groupings of lines are for 
different heating depths, and within each grouping of lines, 
there are three calculations: MESA using grey_irradiated 
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surface boundary condition (solid red line), MESA using the 
Fi,-!,-!, surface heating profile (dashed blue line), and CEPAM 
(Guillot & Morel 1995) using the same grey irradiated bound- 
ary condition (dotted black line; kindly provided by Tristan 
Guillot). The lower curves , corresponding to shallow heat- 
ing, use fiducial values (Kih,Ky) - (10^^, 6x 10""') cm^ g"' and 
give a model radius significantly smaller than the observed 
radius. The upper curves , corresponding to deep heating, 
use (ATthj/Cv) = (10^^,6 X 10"^)cm^g"', yielding significantly 
hotter temperatures deep in the surface radiative zone, which 
slow the cooling enough to agree with the observed radius. 
The choice = 2/Ky gives agreement between the grey ir- 
radiated and Fi,-'Zi, methods, where the factor of 2 accounts 
for the fact that the grey irradiated boundary condition has 
some heating below 2 = l/zcy The radii are at the rth - 2/3 
photosphere for a vertical path into the atmosphere. 

The agreement between all three methods is excellent, at 
the 1-2% level after 100 Myr. The remaining discrepancy 
between the MESA and CEPAM grey irradiated results are 
likely due to different opacity tables, with the MESA result 
using an update of Freedman et al. (2008) (Freedman 2011, 
priv. comm.) while the CEPAM run uses the Allard et al. 
(2001) COND table. The differences at ages < 100 Myr are 
due to different starting conditions. The CEPAM calculation 
started with initial radius 2Rj, whereas the MESA calcula- 
tions started with 5 Rj. The MESA grey irradiated and 
calculations differ at < 100 Myr, likely because the former has 
a fixed thermal opacity while the latter allows the opacity to 
change. 

2.5. Low-Mass Main Sequence Stars 

Most of MESA star's capability to evolve low-mass (M < 
2 Mq) stars was demonstrated in Section 7. 1 of Paper I. MESA 
has seen use in the asteroseismology of helium core flashing 
stars (Bildsten et al. 2012) and the discovery of a new instabil- 
ity from the onset of ^He burning (van Saders & Pinsonneault 
2012). We expect the future use of MESA star for asteroseis- 
mic investigations of these stars to be substantial (see §3). 

The derivation of accurate planetary radii based on tran- 
sits requires accurate radii of the host stars; this motivates 
MESA star investigations of low-mass stars (Lloyd 2011). 
Figure 4 shows 1 and 5 Gyr isochrones at solar composi- 
tion (Y = 0.27, Z = 0.019) from MESA star (solid hnes) 
and Dotter et al. (2008, dashed lines) in the mass-radius di- 
agram. Data points shown in Figure 4 are taken from Torres 
et al. (2010), Caiter et al. (201 1), Irwin et al. (201 1), and Bass 
et al. (2012). This figure is a reproduction of the upper panel 
of Figure 11 from Bass et al. (2012). Figure 4 indicates that 
MESA star is capable of producing mass-radius relations for 
main sequence stars that are consistent with other widely-used 
models as well as observational data. The MESA star models 
were computed using, as much as possible, the same physical 
assumptions as the models used by Dotter et al. (2008). The 
main difference is the equation of state, for which Dotter et al. 
(2008) used FreeEOS'" and MESA star uses a combination 
of the OPAL (Rogers & Nayfonov 2002) and SCVH EOS for 
thermodynamic parameters relevant to this diagram. 

3. ASTEROSEISMOLOGY 

With its highly configurable output options, and its abil- 
ity to calculate asteroseismic variables, MESA star can read- 
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Figure 4. Stellar isochrones at solar composition spanning 0.1 to 1 Mq 
from MESA star (solid lines) and Dotter et al. (2008, dashed lines) in the 
mass-radius plane. The data points plotted are the same as shown by Bass 
et al. (2012). 



ily produce models suitable for use with a range of oscilla- 
tion codes. In addition to its own text output files, MESA can 
produce outputs in formats widely used by stellar oscillation 
codes, such as fgong and osc (Monteiro 2009). 

In Figure 5 we show the evolution of a 1 Mq model 
in the Hertzsprung-Russell Diagram (HRD) and in Tc-pc 
space. These were evolved following the test case found in 
lM_prejns_to_wd, which was modified to include diffusion. 
This runs without user intervention from pre-main sequence 
to white dwarf. To demonstrate the changing stellar struc- 
ture as the model evolves from the main sequence to post 
helium-core burning on the Asymptotic Giant Branch (AGB), 
we show in Figure 6 some of the fundamental quantities ex- 
tracted from the corresponding profile. data files for the 
models marked in Figure 5. These include the Lamb and 
Brunt- Vaisala frequencies defined respectively by 



r 



1 dlnP 



dlnp 



Fi dlnr dlnr 



(3) 
(4) 



where Cj is the adiabatic sound speed and ( is the spherical 
harmonic degree. 

3.1. The Solar Sound Speed Profile 

The seismic properties of the Sun provide a test of stel- 
lar evolution models, and an opportunity to calibrate (Tmlt 
for any particular set of input physics and other assump- 
tions. The MESA star test case solar _calibration pro- 
duces a calibrated Standard Solar Model. Figure 7 shows the 
difference between the helioseismically-inferred solar sound 
speed profile and this model. We also show "Model S" from 
Christensen-Dalsgaard et al. (1996). Both models employ 
comparable input physics and assume solar abundances from 
Grevesse & Noels (1993) and Grevesse & Sauval (1998). One 
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Figure 5. Hertzsprung-Russell diagram and T^-pc evolution of a 1 M© model evolving from pre-main sequence to the white dwarf cooling sequence. The number 
labels denote selected models, for which we show internal profiles in Fig. 6. 



clear improvement since Paper I is a smoother sound speed 
profile at small r/R, which is primarily due to improvements 
in the diffusion module. This is particularly important for 
asteroseismology, where sharp features in the sound speed 
profile can influence the stellar oscillation frequencies. The 
results are based on the solar calibration test case compiled 
with the GNU Fortran compiler version 4.7.2 on Mac OS X 
10.7.5; Appendix B.9 provides information about how the so- 
lar calibration results may depend on different operating sys- 
tems and compilers. 

3.2. New Asteroseismic Capabilities in MESA 

The "astero" extension to MESA star implements an in- 
tegrated approach that passes results automatically between 
MESA star and the new MESA module based on the adiabatic 
code ADIPLS (Christensen-Dalsgaard 2008a, June 201 1 re- 
lease). The MESA module ADIPLS also supports independent 
use for post-processing. 

This astero extension enables calculation of selected pul- 
sation frequencies by MESA star during the evolution of the 
model. This allows fitting to the observations that can include 
spectroscopic constraints (e.g., [Fe/H] and Teff), asteroseismic 
constraints, such as the large frequency separation (Av) and 
the frequency of maximum power (v^ax), and even individual 
frequencies. A variety of approaches for finding a best-fitting 
model are available, including grid searches and automatic 
minimization by the Hooke-Jeeves algorithm (Hooke & 
Jeeves 1961) or by the "Bound Optimization BY Quadratic 
Approximation" (BOBYQA; Powell 2009) technique. These 
searches are user controlled through a number of parameter 
bounds and step sizes. Users also have full control over the 
relative weight assigned to the seismic and spectroscopic parts 
of the;i^^ statistic. 

For the automated x~ minimization, astero will evolve a 
pre-main sequence model from a user defined starting point, 
and find the best match along that single evolutionary track. 
The code then recalculates the track, again initiated at the pre- 
main sequence, with different initial parameters such as mass, 
composition, mixing length parameter and overshoot, and re- 
peats until the lowest has been found. 

Calculating specific mode frequencies is computationally 
intensive. Hence, a number of options exist to improve the 
efliciency of the minimization when individual frequencies 



are included. Bounds can be established on stellar parame- 
ters (e.g., Teff, central H mass fraction, Av), so that ADIPLS is 
invoked only when the model falls within these bounds. This 
enables certain evolutionary stages to be skipped when other 
observational diagnostics rule them out — if a star is known 
to be a red giant, for instance, there is no sense in invok- 
ing ADIPLS when models are on the main sequence. The 
large frequency separation, Av, of the model is calculated 
as the inverse of the sound travel time through the star, Av 
= [2 / dr/c,]-^ (Tassoul 1980; Gough 1986). There is also the 
option to derive Av using simple solar scaling: Av oc (M/R^Y*-^ 
(Kjeldsen & Bedding 1995). 

Moreover, hierarchical approaches to the frequency fitting 
can be selected, saving large amounts of computational time. 
In one case the radial modes are first calculated, and only 
when they match reasonably well are the non-radial mode 
frequencies derived and included in the x^- This is partic- 
ularly beneficial for red giants where the calculation of the 
non-radial frequencies is extremely time consuming. Another 
example is when the time steps in the stellar evolution calcula- 
tions are too large to find an accurate minimum of;^'^. Hence, 
as a further option to increase efficiency while attaining accu- 
racy, the time steps can be set to automatically reduce when 
the model comes close to the "target box" of the observational 
constraints. As for other modules used in MESA star, astero 
offers a range of graphical outputs including an echelle dia- 
gram where the fitting process can be followed in real time. 

There is also an option for including corrections to the 
model frequencies on-the-fly to compensate for the inade- 
quate modelling of the near surface layers of the star. The 
effect, known as the "surface term," is seen as a frequency 
dependent offset between the modelled and observed acoustic 
frequencies of the Sun (e.g. Christensen-Dalsgaard & Thomp- 
son 1997). The offset increases towards higher frequencies 
and is well described by a power law (Kjeldsen et al. 2008). 
MESA star follows the approach described by Kjeldsen et al. 
(2008) for correcting the surface term. 

To illustrate the performance of astero, we show here a fit 
to the star HD43985. The input frequencies and the spectro- 
scopic constraints are from Deheuvels et al. (2010). We first 
ran a wide-range grid search over M, ffMLT. [Fe/H], and Y, in- 
cluding only [Fe/H], T^ff, and Av as observational constraints. 
The results of this initial search guided our starting parame- 
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Figure 6. Internal structure of the six points (indicated by tlie numbers in each panel) marked in Figure 5. The top panel for each points shows luminosity, 
temperature, mass, nuclear energy generation rate, and hydrogen and helium mass fractions. Grey areas mark convective regions according to the Schwarzschild 
criterion. Bottom panels show N and S( for harmonic degrees ( = 1.2, and 3. The dashed line indicates the frequency of maximum power, v„,ax of the 
stochastically excited solar-like modes. 



ters and ranges for the next automatic minimization. We 
first compare our grid results with those of the RADIUS grid 
search routine (Stello et al. 2009), which is based on a grid 
of ASTEC models (Christensen-Dalsgaard 2008b) and find 
agreement within uncertainties. 

We then include the individual oscillation frequencies and 
use the Hooke-Jeeves algorithm for the minimization. 
Model frequencies were corrected for the surface term, and 
the part of the coming from the frequencies was given 2/3 
of the weight in the final similar to that used by Metcalfe 



et al. (2012). To ensure we adequately sample the parameter 
space, we initiate the search at several initial values within a 
broad range. 

Each "Hooke" search generates several stellar evolution 
tracks, each with a best x value. We then combine the data 
from about 1400 tracks to estimate the 1-cr uncertainties in the 
varied parameters following the approach by Deheuvels et al. 
(2010). The lowest (reduced) x^ value we obtained was 2.4 
with a few tens of models in the 2.4—4.0 range, which all fit 
the frequencies similarly well. Among these models there are 
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Figure 7. Comparison of the difference between the helioseismically- 
inferred sound speed profile (Bahcall et al. 1998) to that a MESA star model 
and Model S (Christensen-Dalsgaard et al. 1996). 



two families of results, one of which has slightly lower [Fe/H] 
and Y, and a slightly increased value for the spectroscopic part 
of the x^- 
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HD43985. Observed frequencies are shown with filled symbols as blue 
squares {( = 2), black circles (C = 0), and red triangles (t = 1), and the 
matched model frequencies are shown with open symbols. Black horizontal 
lines indicates l-o" error bars. 

The comparison of the observed and modeled frequencies 
for the realization with the lowest is shown in the echelle 
diagram format in Figure 8. A plot of the internal structure 
including the Brunt-Vaisala and Lamb frequencies is shown 
in Figure 9, and the parameters of the model are listed in Ta- 
ble 2. We set foy - 0.015 and use the GN98 solar abundances. 
Our results can be best compared to those listed as "low aov" 
and "GN93" in Table 4 of Deheuvels & Michel (2011) and 
agree within the uncertainties. 




Figure 9. Same format as in Figure 6, but for the best-fitting model of 
HD49385 (see also Table 2). 



Table 2 

Properties of best fitting 
model to HD49385 



Quantity 



Value 



M/Mo 


1.30 ±0.04 


R/Re 


1.972 ±0.016 


L/Lo 


4.9 ± 0.4 


logg 


3.962 ± 0.003 


Teff/K 


6115 ± 125 


Age/Gyr 


4.1 ±0.4 


^MLT 


1.9 ±0.1 


[Fe/H], 


0.15 ±0.04 


[Fe/H]/ 


0.063 


^initial 


0.29 ± 0.02 


^initial 


0.0222 




2.40 



" [Fe/H]j is the log of the ra- 
tio of the surface (Z/X) rel- 
ative to the solar value of 
0.02293. 

3.3. The Effect of Composition Gradients on the 
Brunt-Vdisdld Frequency 

For g-modes, the most important quantity is the Brunt- 
Vaisala or buoyancy frequency A^. It is crucial for pulsation 
studies that a smooth and accurate method is available for cal- 
culating A^^. In a highly degenerate environment, the pressure 
is nearly independent of temperature and P oc so from 
eq. (4) we see that A^- depends on the difference of two large 
and nearly equal quantities. This can lead to a loss of pre- 
cision and a noisy A^^. To eliminate this problem, A^^ is re- 
written into a form that depends on the difference of the adi- 
abatic and true temperature gradients and on the composition 
gradient: 



A^- 



2 _ 

P Xp 



(Vad -Vt+B). 



(5) 



The term B explicitly takes into account the effect of com- 
position gradients and is commonly called the Ledoux term 
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(e.g., Unno et al. 1989; Brassard et al. 1991). For the general 
case of an A^-component plasma with mass fractions {Z,}, the 
standard formula for B is (e.g., Unno et al. 1989) 



B 



1 '^-^ 

■^T ,= 1 



dlnP 
dlnX: 



P,T.[Xj, 



dlnXi 
dlnP' 



(6) 



Since X, = 1, one of the mass fractions can be elimi- 
nated, so that the sum in eq. (6) runs from 1 to - 1 . We note 
that the partial derivatives in eq. (6) hold all the {Xj} constant 
except for X, and X^, where X^ is varied so as to maintain 

Although eq. (6) is formally correct, we have developed a 
new prescription that is both numerically robust and simple to 
implement. We define a new Ledoux term by taking a varia- 
tion along the radial composition gradient in the model space, 

B- ^ lim l"^0D>7'>^ + (d^/dlnf)51nP)-lnP(p,r,X) 

(7) 

In a stellar model, the implementation of the derivatives in- 
volve differencing quantities on neighboring mesh points. 
Thus, using the subscript k to denote the value of a given 
quantity on the kth mesh point, we have 



B 



1 lnP(p,,n,X,^i)-inP(pk,n,Xk) 



lnP,+i -InP, 



(8) 



This is the form of the Ledoux term that is implemented in 
MESA and we term it the "New Ledoux" formulation. Since 
MESA ensures that Y!i=\ = 1 at each mesh point, this condi- 
tion does not have to be separately enforced. This formulation 

requires just one numerical difference along X that is consis- 
tent with the stellar model and equation of state and does not 
require derivatives with respect to individual mass fractions. 
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Figure 10. A comparison of the new Ledoux prescription for versus the 
direct numerical calculation. This calculation is for a 0.535 Mq white dwarf 
model at T^g = 12, 300 K. 

In Figure 10 we compare A^^ obtained with the New Ledoux 
prescription (grey heavy curve) to that found from the direct 
numerical evaluation of eq. (4) using finite differences (thin 
black curve). For the 0.535 white dwarf model presented 



here it is evident that the new Ledoux criterion gives notice- 
ably smoother results. In particular, the inset box highlights 
a region in which the direct calculation exhibits many noise 
that the new formulation does not. 

4. MIXING MECHANISMS INVOLVING 
COMPOSITION GRADIENTS 

We described the implementation of mixing-length theory 
(MLT) in Paper I, including the allowance for overshoot be- 
yond the boundaries of the convective zones as determined by 
the standard Schwarzschild condition, Vrad > ^ad- Overshoot- 
ing is implemented via an exponential decay of the convective 
diffusion coefficient beyond the boundary of full convection, 
following Herwig (2000): 



Dow - £*conv,o exp 



2Ar 

/ov 



(9) 



where Dcom,a is the diffusion coefficient at the convective bor- 
der, Ar is the distance from the start of overshoot, and Ap is the 
local pressure scale height. The user-adjusted parameter /ov 
then determines the extent of the overshooting region. MESA 
also allows for the adoption of a step-function overshooting 
model, where the mixing region extends a distance /ov/lp be- 
yond the convective boundary with a constant specified diffu- 
sion coefficient. 

In Paper I we did not implement the influence of compo- 
sition gradients on mixing and the resulting diffusion coef- 
ficients when instabilities are operative. The description of 
how MESA star calculates the Ledoux criterion is in §3.3. In 
this section, we describe the implementation of mixing due to 
composition gradients in stellar interiors. 

4.1. Semiconvection 

Semiconvection refers to mixing in regions unstable to 
Schwarzschild but stable to Ledoux, that is 



Vad < Vr < Vl 



(10) 



where Vl is the sum of the adiabatic gradient and the Brunt 
composition gradient term (see eq. [6]), 



Vl = '^,i + B 



(11) 



Once Vl is calculated, regions satisfying equation (10) un- 
dergo mixing via a time-dependent diffusive process with a 
diffusion coefficient calculated by the mlt module following 
Langer et al. (1983), 



Dsc - <3's( 



K 



Vr - Va, 



6Cpp j Vl - Vj- ' 



(12) 



where K = 4acT^ /(3kp) is the radiative conductivity, Cp is 
the specific heat at constant pressure, and asc an efficiency 
parameter 

We stress that semiconvection and overshooting have dis- 
tinct implementations in MESA. Both are time-dependent dif- 
fusive processes. As an example, in Figure 1 1 we display pro- 
files of thermodynamic gradients and their resulting diffusion 
coefficients during core helium burning in a semiconvective 
model with = 0.02 and in an exponentially overshooting 
model with /ov = 10"^. 

4.2. Thermohaline Mixing 
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Figure 11. Sample profiles of semiconvective (left) and exponentially over- 
shooting (right) 3Mo models undergoing core helium burning. Top panels 
show the radiative, adiabatic, temperature, and Ledoux gradients that deter- 
mine mixing boundaries and diffusion coefficients. Bottom panels show the 
resulting diffusion coefficients for energy and cheinical transport. In either 
case, a thin dotted line spanning a single intermediate cell joins the convec- 
tive and semiconvective/overshoot curves. This is intended merely as a guide 
for the eye, as dift'usion coefficients are defined only at the two boundaries of 
a cell. In particular, diffusion for this intermediate cell is governed by con- 
vection at its interior boundary and semiconvection/overshoot at the exterior 
The seiTuconvective model shown here was run with o-sc = 0.01; the expo- 
nentially overshooting inodel with /ov = 10"^. The profiles are taken at the 
points marked in Figure 15. 

Thermohaline mixing arises in the presence of an inversion 
of the mean molecular weight in regions that are formally sta- 
ble against convection according to the Ledoux criterion. 



V„rf < B < 0, 



(13) 



In MESA thermohaline mixing is treated in a diffusion ap- 
proximation, with a diffusion coefficient motivated by the lin- 
ear stability analysis of Ulrich (1972) and Kippenhahn et al. 
(1980) 

3K B 
Ah = ffth (14) 

IpCp (Vt- - Vad) 

The quantity o-th is an efficiency parameter. In the linear anal- 
ysis it depends on the aspect ratio of the blobs/fingers arising 
from the instability. In the case of salt fingers such a value 
is calibrated using laboratory experiments in water (e.g. Kr- 
ishnamurti 2003), where the fingers have an aspect ratio of 
^ 5. In the stellar case the value of this parameter is vexatious 
(e.g. Charbonnel & Zahn 2007; Denissenkov & Pinsonneault 
2008; Cantiello & Langer 2010; WachHn et al. 201 1), with re- 
cent 2D and 3D hydrodynamical calculations pointing toward 
a much reduced value of Cth relative to the salt fingers case 
(Denissenkov 2010; Traxler et al. 2011). Figure 12 shows a 
calculation including the effects of thermohaline mixing dur- 
ing the RGB phase of a IM© star after the luminosity bump 
(e.g. Chai-bonnel & Zahn 2007; Cantiello & Langer 2010). 

4.3. Impact of Mixing on Convective Core Hydrogen and 
Helium Burning 



The duration of the hydrogen and helium core burning de- 
pends on the extent of the convective core, so we focus here 
on exhibiting the MESA capabilities during these phases. As 
we noted above, there are many physical effects that change 
the size of the convective core, such as semiconvection, over- 
shooting, and rotation-induced mixing. For example, the 
Schwarzschild criterion implies larger cores than the Ledoux 
criterion, but when using Ledoux alone, the region above the 
convective boundary is overstable and so semiconvection oc- 
curs (see §4.1). 

We evolved a non-rotating 1 .5 Mg star with 
(y, Z) = (0.23,0.02) through central hydrogen burning 
using Ledoux, Ledoux plus semiconvection, Schwarzschild, 
and Schwarzschild plus overshoot. As is evident in Figure 
13, this set of physical processes leads to a large range of 
convective core masses and thereby main sequence lifetimes. 
For the parameters explored we found that overshooting 
increases the lifetime by a factor <1.2 for Schwarzschild 
and <2.5 for Ledoux. Figure 14 shows an HR diagram 
for each of the 1.5 M© models undergoing core hydrogen 
burning, showing the impact of convective core extent on 
main-sequence turnoff morphology. 

We also evolved a non-rotating 3Mq star with {Y,Z) - 
(0.25, 0.02) through central helium burning. Overshooting ex- 
tends the burning lifetime by a factor <1.6 for Schwarzschild 
and <2.8 for Ledoux. Although this lengthening of the core 
burning phase is always true of convective overshoot, we find 
that the extension of the overshoot and convective regions is 
sensitive to the temporal resolution adopted. With sufficiently 
large values of /ov the upper boundary develops oscillatory 
behavior which can also affect the lifetime. This behavior 
also occurs when the overshoot is implemented with the step- 
function implementation of overshoot. This instability is not 
seen in overshoot during hydrogen burning and has yet to be 
studied in detail. 

5. EVOLUTION BEYOND THE MAIN SEQUENCE AND 
WHITE DWARFS 

Extending the verification of Paper I, we now compare to 
other available codes for intermediate-mass stars, 3-8 M©. 
We describe the techniques used by MESA star to evolve stars 
through the AGB phase to the white dwarf cooling sequence. 
We also demonstrate how MESA star incorporates compres- 
sional heating from accretion. 

5.1. Code Comparisons during Helium Core Burning 

We start by comparing the results of MESA star to those 
from the Dartmouth Stellar Evolution Program (DSEP; Dot- 
ter et al. 2008) for stars with M - 3-8 Mg. In both cases, 
the models were evolved from the pre-main sequence to the 
depletion of helium in their cores. For completeness, the 
MESA star models were further evolved to the occurrence of 
the first helium thermal pulse. 

All models have an initial composition Y - 0.272, Z - 
0.02, and no mass loss or rotation was included. The bound- 
aries of mixing zones are determined by the Schwarzschild 
criterion with ckmlt - 2. In order to compare the codes, we 
do not allow overshooting or semiconvection. We adopt the 
Kunz et al. (2002) rate for '^C(q', 7)'*0 and the Imbriani et al. 
(2004) rate for ''*N(/5,y)'^0; for all other rates we use the 
NACRE compilation (Angulo et al. 1999). We use the OPAL 
Type 2 opacity tables (Iglesias & Rogers 1993) to account 
for the carbon- and oxygen-enhanced opacities during helium 
burning. 
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Figure 12. Thermohaline mixing during the RGB phase of a Z = 0.02, 1 M© model, initially rotating with an equatorial velocity of 10 km . In the left panel 
a Kippenhahn diagram shows, in mass coordinate and as function of model number, the locations of the retreating convective envelope (blue), of the H-burning 
shell (red) and of the thermohaline mixing region (magenta). The right panel shows diffusion coefficient profiles extracted at model number 1849, which is 
the last model shown in the Kippenhahn plot. The H-burning shell and the convective envelope are shaded in red and blue, respectively. Thermohaline mixing 
(magenta line) transports chemicals between the burning shell and the convective envelope. Also shown are the diffusion coefficients resulting from Eddington- 
Sweet circulation (ES), magnetic torques by dynamo generated fields (ST), Dynamical Shear (DSI), Secular Shear (SSI) and Goldreich-Schubert-Fricke (GSF) 
instability (see §6 for details). 
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Figure 13. History of convective core extent during the main sequence for 
a non-rotating I.SMq star with various mixing options. The plot shows the 
boundary of convection not including the extent of semiconvection or over- 
shooting. 



The resulting tracks in the HR diagram of Figure 16 and 
the evolution in the T^-pc plane of Figure 17 show excellent 
agreement between the codes. Figures 18 and 19 show the 
hydrogen-burning luminosity, the helium-burning luminosity, 
and the extent of the convective core during convective helium 
core burning for a 4Mo model (Fig. 18) and a 6Mq model 
(Fig. 19). Table 3 gives a summary of the core hydrogen burn- 
ing lifetime, the core helium burning lifetime, the final extent 
of the convective core during central helium burning, and the 
final carbon mass fraction Xq in the core for each model. For 
the MESA models, we also show the maximum extent of the 
convective core during central hydrogen burning, the mass of 
the helium core before helium ignition, and the mass of the 
C/O core at the time of the first hehum thermal pulse. 



Figure 14. The HR diagram for the non-rotating 1 Mq star with various mix- 
ing options. Tracks are displayed from ZAMS until depletion of core hydro- 
gen to X = 10"^ 

We close with an additional comparison of the helium core 
burning phase of a M = 3Mq, Z = 0.02 model computed 
by MESA to that of Straniero et al. (2003). Both models were 
evolved using the Kunz et al. (2002) rate for ^^C(a,yy^O. 
The results for MESA star are a helium core burning lifetime 
of 83.6 Myr and final C/O mass fractions of Xq - 0.43, Xq - 
0.55; Straniero et al. (2003) find a lifetime of 88 Myr and 
Xc = 0.42, Xo = 0.56. 

5.2. Making and Cooling White Dwarfs 

In the previous section, we discussed the evolution of 
3-8 Mq stars up to the occurrence of the first He thermal 
pulse. In Paper I we showed detailed comparisons of the 
evolution of a 2 Mq star to the EVOL code (Herwig 2004), 
exhibiting the ability of MESA star to calculate multiple he- 



Modules for Experiments in Stellar Astrophysics (MESA) 



13 



Properties of the 3- 
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! Mq evolution (masses in solar units). Selected quantities are also shown from DSEP for comparison. 
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Figure 15. History of convective core extent during the core helium burn- 
ing phase for a non-rotating 3Mo star with various mixing options, as in 
Figure 13. Time is measured relative to the onset of the convective core 
burning. Efficient semiconvection (age = 0.01) and inefficient overshooting 
(fay = 10"^) coincide with the pure Schwarzschild model. The filled (open) 
circle indicates the time for which we display a profile detailing semiconvec- 
tion (overshooting) in Figure 1 1 . 



Hum shell pulses. We now illustrate the final evolution of 
intermediate-mass stars, and how to construct white dwarfs 
(WDs) by using winds. 

We evolve 3 Mq, 5 Mq, and 7 Mq stars from the ZAMS us- 
ing the test suite case make_co_wd. This makes use of RGB 
mass loss following Reimers (1975) with an efficiency param- 
eter T] - 0.5 and AGB mass loss following Bloecker (1995) 
using = 0.1 until the occurrence of the first helium shell 
flash. At that time, an increased Bloecker 77 = 5 is adopted to 
allow only a small number of thermal pulses before the wind 
mass loss eliminates the envelope. Such intervention allows 
MESA star to make a high-mass WD. To avoid shortening of 
timesteps due to radiation-dominated envelopes, these cases 
also use the MLT+-I- capability described in §7.2. 

Figure 20 shows the resulting tracks on the HR diagram. 
The 3 Mq star underwent eight thermal pulses after the en- 
hancement of Bloecker winds, while the 5 Mq and 7 Mq stars 
lost their envelopes so quickly that thermal pulses were imme- 
diately halted. The 5Mo star ended up as an M = 0.844 Mq 
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Figure 16. Hertzsprung-Russell diagram for evolution of 3-8 Mq stars from 
the pre-main sequence through core helium depletion. Models are from MESA 
(thick grey lines) and DSEP (dashed black lines). Each curve is labeled with 
its corresponding initial mass in solar units. 

C/0 WD with a helium shell of thickness Mee = 0.009 Mq 
and a hydrogen envelope of Mh = 2.3 x 10"^ Mg. Note that 
the C/O WD mass is only slightly larger than the C/O mass at 
the first thermal pulse (0.827 Mq) reported in Table 3. 

After removal of the envelope, the evolution of the white 
dwarf is continued through its cooling phase past solidifica- 
tion. We include gravitational settling and chemical diffu- 
sion of the outermost layers. Figure 21 shows T-p profiles 
taken at various effective temperatures during the cooling of 
the M = 0.844 Mq C/0 WD made from the 5M0 star The 
growing depth of the convection zone is shown by the dashed 
line, and the open circles designate the H/He transition, while 
the filled circles denote the He/CO transition. Figure 22 illus- 
trates the resulting L-T^ relation as these models cool. 

The test suite case wd_di£fusion uses the implementation 
of diffusion described in Paper 1 to evolve a WD of mass 
0.611 Mq until the Me = 1.3 x IO^'^Mq hydrogen layer and 
the Mhe - 4.6 X 10"^ Mq helium layer approach diffusive 
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Figure 17. Same as Fig. 16, but in tlie T^-pc plane. The MESA models (thick 
grey lines) are evolved until the occurrence of the first thermal pulse. 
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Figure 18. History of hydrogen burning luminosity (top), helium-burning lu- 
minosity (center), and convective core extent (bottom) during the core helium 
burning phase for the 4 M© models. Time is measured relative to the onset of 
the convective core. 



equilibrium. At this point, the WD has an effective temper- 
ature of Teff ~ 5,000K. We show the resulting abundance 
profiles in Figure 23, and, for comparison, the abundance pro- 
files derived from the analytic form for diffusive equilibrium 
(eq. (22) of Althaus et al. 2003). This formula is obtained by 
integrating equation (A. 5) of Arcoragi & Fontaine (1980) and 
assuming an ideal gas equation of state and complete ioniza- 
tion of both species. 

The specific treatment of convection can also impact WD 
evolution. In Paper 1, MESA used the Cox & GiuU (1968) pre- 
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Figure 19. History of hydrogen burning luminosity (top), helium-burning lu- 
minosity (center), and convective core extent (bottom) during the core helium 
burning phase for the 6 Mq models. Time is measured relative to the onset of 
the convective core. 



L I I I I I I I I I I I I I I I I I I I I I I 





I— 1 — 

O 



(N 

I 



a. - 




"lllllllllllllllll^llllll" 



5.5 5.0 4.5 4.0 

logTeff (K) 



3.5 



Figure 20. Evolution of 3. 5 and 7 Mq models from zero-age main sequence 
to cooling white dwarfs. A Bloecker mass loss strips the stars of their en- 
velopes on the thermally pulsing AGB to make the three C/O white dwarfs. 
The single 0.32 Mq He white dwarf was made with mass loss after the hydro- 
gen main sequence for the 3 Mq model was completed. 



scription for convection as its default convective MLT, with 
the optional extension of Henyey et al. (1965). Since Paper 1, 
we have added support for the formulations of Bohm-Vitense 
(1958), Bohm & Cassinelli (1971), and Mihalas (1978). In 
particular, the Bohm & Cassinelli prescription, often referred 
to as "ML2," is frequently employed in WD studies (e.g., 
Bergeron et al. 1995). In Figure 24 we show a comparison of 
the Brunt-Vaisala frequency calculated with MESA to that us- 
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Figure 21. Profiles in log T-logp space of the cooling 0.844 M© C/O white 
dwarf evolved from a 5 Mq progenitor. Each model is labeled to the right 
by Tiff. The outermost point of the model is at r = 25. Dotted curves de- 
note convective regions. Going toward the interior, open circles designate the 
transition into the helium-rich shell, and filled circles designate the transition 
into the C/O core. 
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Figure 22. Surface luminosity as a function of central temperature for the 
cooling 0.32, 0.574, 0.844, and 0.928 Mq WDs evolved from 3, 5, and 7Mo 
progenitors. 



ing the Warsaw envelope code (Paczyriski 1969, 1970; Pamy- 
atnykh 1999), assuming the ML2 prescription. This is the 
same WD as in Figure 23, but now at a lower Teff = 11, 354 K. 
To more accurately integrate these opaque but thin layers, we 
reduce r at the boundary of the model by a factor of 1000 from 
its photospheric value of 2/3. This calculation is a sensitive 
test of the envelope integrations because A^- is a derivative of 
the envelope structure. The two codes give indistinguishable 
results for this case and all other cases that we have calculated. 

MESA now includes atmospheric tables based on the non- 
grey model atmospheres for hydrogen-atmosphere WDs 
(Rohrmann 2001; Rohrmann et al. 2012). Such an approach is 
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Figure 23. A comparison of time-dependent dilfusion calculations for a 
M = 0.611 Mo WD with M„ = 1.3 X IQ-^Mq and Mhc = 4.6 X IO^-Mq 
with MESA star (solid lines) to those assuming diffusive equilibrium and an 
ideal gas equation of state (dashed lines). 
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Figure 24. A comparison of the Brunt- Vaisala frequency calculated with 
MESA (solid grey line) to that using the Warsaw envelope code pig35.£ 
(dashed line) for the same WD in Figure 23, but at a cooler %fi = 1 1, 354 K. 



necessary at T^ff < 6000 K, where WDs develop deeper con- 
vection zones. When the convection zone comes in contact 
with the degenerate, nearly isothermal core, energy is able 
to flow out of the core much more efficiently. The use of 
non-grey atmosphere models results in shallower convection 
zones, so this convective coupling of the core and envelope 
is delayed. For reliable cooling ages, we therefore recom- 
mend using non-grey atmospheres when Teff < 6000 K. Fig- 
ure 25 demonstrates the impact of non-grey atmospheres with 
the 0.535 M© WD, which was cooled with and without the 
non-grey atmosphere. 

MESA currently treats crystallization by employing the 
Potekhin & Chabrier (2010) EOS (PC EOS). The phase tran- 
sition is first-order, so the EOS exhibits a latent heat between 
the solid and liquid phases, i.e., the entropy and internal en- 
ergy both experience finite jumps. This energy is included in 
MESA star models of cooling white dwarfs through the grav- 
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Figure 25. The age difference (non-grey minus grey) in Gyr as a function of 



itational source term in the energy equation. 



Sera 



dS 
dt 



(15) 



This form for Cg^-^y replaces the default one (see eq. (16) be- 
low) in cells where F > 160, and is smoothly interpolated 
with the default form in cells where 130 < F < 160. The 
PC EOS uses the criterion F - 175 to determine crystalliza- 
tion, but it is straightforward to include explicit crystallization 
curves for C/0 and other mixtures (e.g., Schneider et al. 2012; 
Medin & Gumming 2010). For example, using the parame- 
ters of the model in Figure 23, the age difference at late times 
(Teff < 3,500K) between a model with and without the la- 
tent heat of crystallization is x 0.8 Gyr; a slightly larger value 
would be obtained using the phase diagram of Schneider et al. 
(2012). MESA does not currently treat phase separation of dif- 
ferent chemical species upon crystallization. 

Low mass WDs (M < 0.4 M©) with helium cores and 
hydrogen envelopes may be produced in binary systems 
when the envelope is stripped by the companion as the pri- 
mary evolves up the giant branch (Iben 1991, and references 
therein). He-core WDs of mass M ^ 0.4-0.5 Mq may also 
be produced through strong RGB winds (D'Cruz et al. 1996), 
although we do not discuss this possibility further here. 

Here we discuss the prescription for stripping the envelope 
used in the test case makeJie_wd. The first step is to evolve 
a star, M - 3.0 M© in this example, from the PMS until a 
He core of the correct size has been made. The remnant to- 
tal mass is determined by the mass interior to where the H 
abundance has dropped below a preset value, for example, 
Xu = 0.1, moving in from the surface. Next, the routine 
relaxjnass is used to remove mass from the model until it 
has the desired remnant mass. After the initial remnant has 
been constructed, diffusion can then be turned on to allow an 
outer H layer to form. After this stage, normal evolution of 
the WD occurs, as shown in Figures 20 and 22. 

5.3. Compressional Heating and Accretion 

As discussed in §6.2 of Paper 1, MESA star calculates egrav 
of eq. (15) in terms of the local thermodynamic variables (T 



and p) used by MESA, 

ftrrav — —CpT 



/, „ xdlnT „ dlnp 
(1 - Vad^rj-^ -^adA-p 



df 



(16) 



MESA star takes the quantities in this equation as provided 
by EOS, and computes the Lagrangian time derivatives to find 
fgrav- MESA star can alternatively work under the assumption 
that P = i|as + ^ad, in which case MESA star treats i|as rather 
than p as its basic variable (see §B.2 for a discussion). In that 
case. 



>^grav 



= -CpT 



i?,,d\dlnr 



dt 



- V, 



Pga.dlnPga 



ad" 



dt 



(17) 

Either formulation can be used deep within the star, as long 
as the location is safely removed from any phase transition. 
Paper 1 described the validation of these formulations. 

A complication arises when egiav needs to be evaluated in 
material that was not present in the previous timestep. Defin- 
ing the envelope mass coordinate AM = M - m, we need to 
resolve the entropy for AM < 6M - M 6t, as the explicit La- 
grangian time derivatives of eqs. (16) and (17) cannot be nu- 
merically evaluated. Since there can be important physics that 
needs to be resolved for these mass shells for AM <s: SM, an 
approximation must be derived that allows for accurate mod- 
eling of the star's outermost layers without having to result to 
a dramatic shortening of St. 

The luminosity L^cc = GMM/R from the accretion shock 
(or boundary layer) goes outwards and does not determine 
the entropy of the material as it becomes part of the hydro- 
statically adjusting star Rather, the entropy of the material at 
AM 6M is determined by the the transport of L (Nomoto & 
Sugimoto 1977; Nomoto 1982; Townsley & Bildsten 2004). 
Consider such an outermost layer, where there are two rele- 
vant timescales, the thermal time, fth = CpT AM I L, and the lo- 
cal accretion time, f^cc = AMjM. In nearly all relevant cases, 
the ratio t^\,jt^^^ - CpTMjL <K 1; this implies that the fluid 
element adjusts its temperature to that needed to transport the 
stellar luminosity from deep within. This simplifies e^^-^ in 
that part of the star (following Townsley & Bildsten 2004) to 



CpTGmM 



(Vad - Vr), 



(18) 



enabling accurate modeling within MESA star of nearly 
all fluid elements that become part of the star during each 
timestep, many of which have envelope mass coordinates 
AM <K M6t. 

We give an explicit example of this thin-shell radiative cal- 
culation of egrav in a C/O white dwarf accreting hydrogen- 
rich material and undergoing classical nova (CN) cycles. We 
present two models accreting at rates of M = 10"" Moyr"' 
and 10"'°Moyr-i. Both cases were evolved from a 0.6 Mq 
starting model with T^ - 10^ K which had undergone a few 
flashes while accreting at M = 10"" Mq yr"'. The accreted 
material has solar-like composition X - 0.70, Y - 0.26, and 
Z = 0.04 where the metal mass fractions are taken from Lod- 
ders (2003). 

Profiles of the envelope during the mass accumulation 
phase between CN outbursts for the two accretion rates are 
displayed in Figures 26 and 27. Each line represents a differ- 
ent time in the accumulation cycle up to the unstable ignition, 
when the hydrogen mass reaches Me - Mign. All material at 
pressures smaller than that shown by the open circle is new to 
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the model in that timestep (e.g., it has AM < 5M) and employs 
the modified egrav of eq. (18). This highlights the significance 
of this approximation as it allows MESA star to calculate ma- 
terial properties at AM ~ 10"** 6M. The solid points show 
where e^^v switches to the explicit form employing the La- 
grangian time derivatives, such as eq. (16). 

The middle panel shows egiav^* fgiavAM, which reflects 
the contribution of egmv to the outward luminosity. The dis- 
continuity of egiav at the solid point reflects the error associ- 
ated with the abrupt transition in the calculational approach. 
The substantially larger luminosity of the early (Me/Mign = 
0.22) stages is due to the ongoing transfer of heat from the pre- 
vious outburst. The near-discontinuous drop in egrav occurs at 
the base of the hydrogen-rich envelope, and reflects the jump 
in composition from the accreted material to the nearly pure 
^He layer The expected amplitude of the jump in egrav de- 
pends on both the composition jump and the local degree of 
electron degeneracy (see Appendix B of Townsley & Bildsten 
2004 for a discussion). 




log P (erg cm ) 

Figure 26. Envelope profiles as a function of pressure of the accreting 
white dwarf for three instants during the mass accumulation phase; M = 
10"" Moyr"' model. The top panel shows temperature, the central panel 
shows the gravitational energy release rate, and the bottom shows the lumi- 
nosity. Material to the right of the open circle is newly accreted. The code 
treats material to the right of the filled circle using the thin-shell radiative 
calculation of e„nv 



6. ROTATION 

A star's rotational energy is usually a small fraction of 
the gravitational energy: for the Sun it is ~ 10"^ and for a 
25 Mq star rotating with a typical equatorial velocity - 
200 km s"' on the main sequence it is ~ 0.04. Therefore 
the effects on the stellar hydrostatic equilibrium are marginal, 
with the exception of stars close to critical rotation (see §6.3). 
Even in the case of a small perturbation to hydrostatic equi- 
librium, rotation induces a modification to the star's thermal 
equilibrium (von Zeipel 1924). Together with the emergence 
of rotationally-induced dynamical and secular instabilities, 
this can significantly affect the evolution of stars (Maeder & 
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Figure 27. Same as Fig. 26, but for model accreting 10"'" Mq yr"' . 

Meynet 2000b). Due to the destabilizing effect of increasing 
radiation pressure, rotation is particularly important in mas- 
sive stars (see, e.g., Heger et al. 2000; Meynet & Maeder 
2000). Moreover, the final fate of a massive star depends 
chiefly on the relative importance of rotation during its evolu- 
tion (e.g., Heger et al. 2000; Hirschi et al. 2004; Heger et al. 
2005; Yoon et al. 2006; Woosley & Heger 2006; Ekstrom et al. 
2012; Georgy et al. 2012; Langer 2012). 

Here we describe the implementation of rotation in 
MESA star. We briefly discuss the modification to the stel- 
lar structure equations and the inclusion of rotationally- and 
magnetically-induced mixing. Magnetic fields generated by 
differential rotation in radiative regions have been imple- 
mented following the work of Spruit (2002) and in the same 
fashion as in Petrovic et al. (2005) and Heger et al. (2005). 
Rotationally enhanced mass loss is also discussed. 

We compare rotating massive-star models calculated with 
MESA star to previous calculations performed with KEPLER 
(Heger et al. 2005). We also directly compare runs from 
MESA star and STERN (Petrovic et al. 2005; Yoon & Langer 
2005; Brott et al. 2011). The purpose of these tests is to 
verify our implementation of rotation, which is derived from 
STERN. We do not compare to codes that have a different 
implementation of rotation (e.g., Hirschi et al. 2004; Ekstrom 
et al. 2012; Georgy et al. 2012; Potter et al. 2012b,a). Al- 
though beyond the scope of this paper, such comparisons are 
critical when coupled to observations of the effects of rota- 
tion in stars (e.g.. Hunter et al. 2007; Evans et al. 2011) and 
asteroseismology (Beck et al. 2012; Mosser et al. 2012). 

6.1. Implementation of Shellular Rotation 

Stellar structure deviates from spherical symmetry in the 
presence of rotation. While the structure is inherently three- 
dimensional, it suffices to solve the stellar structure equa- 
tions in one dimension if the angular velocity, oj, is constant 
over isobars (the so-called shellular approximation; see, e.g., 
Meynet & Maeder 1997). This is expected in the presence of 
strong anisotropic turbulence acting along isobars. In radia- 
tive regions such turbulence is a consequence of differential 
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rotation (Zahn 1992) and efficiently erases gradients along 
isobars and enforces shellular rotation (Meynet & Maeder 
1997). Turbulence in the vertical direction (i.e., perpendicular 
to the isobars) is much weaker due to the stabilizing effect of 
stratification. In MESA star we adopt the shellular approx- 
imation (Meynet & Maeder 1997) and calculate the modifi- 
cation to the stellar equations due to centrifugal acceleration 
(Kippenhahn & Thomas 1970; Endal & Sofia 1976). 

Transport of angular momentum and chemicals due to 
rotationally-induced instabilities is implemented in a diffusion 
approximation (e.g., Endal & Sofia 1978; Pinsonneault et al. 
1989; Heger et al. 2000). This choice has also been adopted 
by other stellar evolution codes (e.g., KEPLER, Heger et al. 
2000; STERN, Yoon & Langer 2005). We stress that this is 
not the only possibility, and other groups have impzlemented 
a diffusion-advection approach (e.g., GENEVA, Eggenberger 
et al. 2008; RoSE, Potter et al. 2012b). The RoSE code can 
switch between the two different implementations. The two 
approaches are equivalent for the transport of chemicals. Po- 
tentially large differences can arise, however, for the trans- 
port of angular momentum. A detailed description of the 
advection-diffusion equation for angular momentum is given 
in Zahn (1992) and Maeder & Zahn (1998). 

In MESA star the turbulent viscosity v is determined 
as the sum of the diffusion coefficients for convec- 
tion, semiconvection, and rotationally-induced instabilities. 
MESA star calculates diffusion coefficients for five different 
rotationally-induced mixing processes: dynamical shear in- 
stability, Solberg-H0iland instability, secular shear instability, 
Eddington-Sweet circulation, and the Goldreich-Schubert- 
Fricke instability. See Heger et al. (2000) for a detailed de- 
scription of the physics of the different instabilities and the 
calculation of the respective diffusion coefficients. These en- 
ter the angular momentum and abundance diffusion equations 
that are solved at each timestep (see §B.6). 

6.2. Magnetic Fields 

Differential rotation in the radiative layers of a star can am- 
plify a seed magnetic field. Such a dynamo process has been 
suggested by Spruit (2002, Spruit-Tayler dynamo); a debate 
on this is still ongoing (Braithwaite 2006; Zahn et al. 2007; 
Denissenkov & Pinsonneault 2007). Observations of the final 
spins of both WDs and neutron stars (Heger et al. 2005; Suijs 
et al. 2008) suggest that angular momentum transport with an 
efficiency similar to the torques provided by the Spruit-Tayler 
dynamo operates. Models that only include angular momen- 
tum transport through rotational instabilities do not produce 
the core-envelope ratio of angular velocity observed through 
the splitting of mixed modes in red giant stars (Eggenberger 
et al. 2012). 

MESA star accounts for transport by magnetic fields of 
angular momentum and chemicals due to the Spruit-Tayler 
dynamo. We refer to Spruit (2002) for a description of the 
physics of the dynamo loop and to Maeder & Meynet (2003), 
Maeder & Meynet (2004) and Heger et al. (2005) for a discus- 
sion of its inclusion in stellar evolution codes. We implement 
the Spruit-Tayler dynamo in MESA star following KEPLER 
(Heger et al. 2005) and STERN (Petrovic et al. 2005). 

6.3. Rotationally-Enhanced Mass Loss 

We include the rotational modification to the mass loss 
rate (Friend & Abbott 1986; Langer 1998; Heger & Langer 
1998; Maeder & Meynet 2000a). Similar to other codes (e.g., 



Heger et al. 2000; Brott et al. 2011; Potter et al. 2012a), in 
MESA star the stellar mass loss is enhanced as the rotation 
rate increases according to the prescription 



M(Q) = M(0) 



1 



1 — Q/Qcrit 



(19) 



where Q. is the value of the surface angular velocity and Qcrit 
is the critical angular velocity at the surface. This last quan- 
tity is defined as Q?:^.^^ = (1 - L/LEdd)GM/R^, where LEdd - 
AncGMjK is calculated as a mass-weighted average in a user- 
specified optical depth range (default value t e [1 - 100]). In 
MESA star the default value for the exponent ^ is 0.43 (Langer 
1998). Other implementations of rotationally enhanced mass 
loss can be found in Maeder & Meynet (2000a) and Georgy 
et al. (2011). 

For stars approaching Q/Qcrit = 1, the mass loss calculated 
using equation (19) diverges. Notice that luminous stars can 
approach this limit without having to rotate very rapidly as 
f^crit — > when L/^Edd — > 1- Following Yoon et al. (2010) we 
limit the mass loss timescale to the thermal timescale of the 
star tkh 



M - min 



M 

M(n),f 



(20) 



where / is an efficiency factor of order unity (default value is 
/ = 0.3). 

6.4. Test Cases: 15 Mq and 25 Mq 

As a first test we initialize a 15Mq model with Z - 0.02 and 
initial equatorial rotational velocity - 200 km s"' and run 
two calculations: 

• 15MAG includes the effects of rotation and Spruit- 
Tayler magnetic fields on both the transport of chem- 
icals and angular momentum. 

• 15R0T includes only the effect of rotation on both the 
transport of chemicals and angular momentum; 

The initial conditions have been calibrated to match as closely 
as possible the KEPLER 15 Mq models (Heger et al. 2005). 
Moreover, we directly compare the MESA star models with 
calculations from STERN (see e.g., Yoon & Langer 2005; 
Yoon et al. 2006). In particular we adopt a value of fc — \ /30 
for the ratio of the turbulent viscosity to the diffusion coeffi- 
cient and a value =0.1 for the sensitivity to /v-gradients 
(see Heger et al. 2000, for a discussion of these calibra- 
tion parameters). The Ledoux criterion is used for the treat- 
ment of convective boundaries, together with semiconvection 
(ffsc = 1). We use ffMLT - 1-6, mass loss as in Yoon et al. 
(2006) with rotational enhancement as described in § 6.3. 

In Fig. 28 we show the evolutionary track and the evolu- 
tion of surface equatorial rotational velocity for the 15MAG 
model. Results of a similar calculation using STERN are 
shown as a dashed curve. The two results are in excellent 
agreement. Small differences in luminosity and lifetimes are 
not unexpected, as we have only matched the physics of ro- 
tation between the two calculations and not other ingredients. 
Values for the diffusion coefficients for rotationally induced 
mixing and magnetic torques during the main sequence of 
15MAG are shown in Fig. 29. The comparison reveals a very 
good agreement. Both stars are kept in solid-body rotation 
during the main sequence by the efficient transport of angular 
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momentum provided by the Eddington-Sweet circulation and 
Spruit-Tayler magnetic fields. 

The amplitude and location of the azimuthal (B^) and ra- 
dial (B,) components of the magnetic fields during different 
phases of the evolution of 15MAG are shown in Fig. 30. As 
expected, these fields are generated only in radiative regions 
of the star and B^ > B,- (Spruit 2002). As the star evolves 
away from the main sequence its structure departs from solid- 
body rotation with the core rotating faster than the envelope. 
During this stage the role of magnetic fields is very important 
in transporting angular momentum from the core to the en- 
velope. The effect can be seen in Fig. 31, which shows the 
evolution of the internal specific angular momentum in mod- 
els 15ROT and 15MAG. The presence of magnetic torques 
results in a dramatic spin-down of the core of 15MAG with 
respect to 15ROT (see also Table 4). These results are in very 
good agreement with the ones obtained by STERN and KE- 
PLER. 

As a second test, we now evolve a 25 M© model (25MAG) 
with the same physics as in 15MAG. Figure 32 directly com- 
pare results with calculations performed with STERN. In 
Fig. 33 we show a detailed comparison of the evolution of the 
internal specific angular momentum profile. We find a very 
good quantitative agreement between MESA star and STERN 
down to He depletion in the core. The timescale for nuclear 
burning decreases substantially after He-burning and becomes 
shorter than the angular momentum transport timescale after 
C depletion. Thus only minor changes in the final angular 
momentum content of the stellar core are expected after this 
stage. Figure 34 shows the full evolution of the specific an- 
gular momentum profile of the MESA star calculation from 
ZAMS to Si exhaustion. 

6.5. Rapidly Rotating Massive Stars 

MESA star can calculate the evolution to core collapse of 
rapidly rotating massive stars. Rotational instabilities can be 
efficient enough to erase the compositional gradients built 
by nuclear burning. In such cases the model never devel- 
ops a compositional stratification and remains almost com- 
pletely mixed throughout its evolution (Maeder 1987b; Yoon 
& Langer 2005; Woosley & Heger 2006). This process leads 
to a bifurcation in the HR-diagram, with stars above a certain 
mass and rotation rate becoming more luminous and hotter. 
The threshold required for this bifurcation depends mostly on 
the initial mass of the star (Yoon et al. 2006). Metallicity also 
plays an important role, as angular momentum is lost through 
line-driven stellar winds, with mass-loss rates depending on 
the metallicity at the stellar surface (Vink et al. 2001). For 
the calculations in this section, we adopt the same mass-loss 
prescription of Yoon et al. (2006). 

Figure 35 shows the evolution of two 16 Mq models at 
metallicity Z - 0.0002 with rotation initialized at the ZAMS. 
One model is rotating very rapidly, with Ueq - 450 km s"' 
(corresponding to Q/Qcrit = 0.55 and J = 3.23 x 10^^ ergs s), 
while the other rotates at Ueq - 280 km s"' (corresponding 
to Q/Qcrit = 0.4 and J = 2.52 x lO^^ergss). The model 
with n/Ocrit = 0.55 avoids the core-envelope structure and 
becomes a compact Wolf-Rayet star The absence of a RSG 
phase eliminates the large magnetic torques from an extended 
envelope. The evolution of the internal profile of specific an- 
gular momentum in the two models clarify this point: the 
model with Q/Qcdt = 0.4 becomes a RSG, and the core 
spins down rapidly. When it reaches core-collapse its struc- 



ture is extended, as implied by the large free-fall timescale 
shown in the left panel of Fig. 35. As a consequence, there 
is not enough angular momentum in its core to build an ac- 
cretion disk around a newly formed compact object. This 
model is expected to produce a Type IIP supernova. On 
the contrary, the model with Q/Qcrit - 0.55 is compact (the 
free-fall timescale is on the order of seconds, right panel of 
Fig. 35) and has enough angular momentum to produce an 
accretion disk around the central compact object. Therefore 
this model is a candidate progenitor for a long gamma-ray 
burst (Woosley 1993). This last calculation can be directly 
compared to the KEPLER model 16T1 in Woosley & Heger 
(2006). 

We further test MESA capabilities by evolving two rotat- 
ing 40 M0 models at Z = 10"^. One model is initialized 
at the ZAMS with - 260 km s"', while the other has 
feq = 630 km s"'. The results of these calculations can be 
compared with the models shown in Yoon & Langer (2005). 
Figure 36 shows that for the more rapidly rotating model, ro- 
tational mixing (mainly due to the Eddington-Sweet circula- 
tion) is large enough that the star evolves blueward in the HR- 
diagram. This evolution results in a compact configuration 
and enough angular momentum to fulfill the requirements of 
the collapsar scenario for long gamma-ray bursts, as shown 
in Fig. 37 (right panel). On the other hand, the slower rotat- 
ing model becomes a RSG and loses most of its core angular 
momentum, as shown in Fig. 37 (left panel). 

7. MASSIVE STELLAR EVOLUTION 

Modeling massive stars is numerically difficult. One prob- 
lem is they develop loosely bound, radiation pressure domi- 
nated envelopes that can cause density and gas pressure inver- 
sions. This environment poses a physical and numerical chal- 
lenge that all stellar evolution codes must address to evolve 
massive stars past the main sequence. In this section we dis- 
cuss MESA star's capability to evolve rotating massive stars 
from their zero age main sequence to core-collapse. 

7.1. Evolution of Massive Stars with MESA 

Previous computations with MESA star found these en- 
velopes to be numerically (and probably physically) unsta- 
ble. This is a known issue in the literature (e.g., Maeder 
1987a), which reveals the limitations of the ID treatment 
of late phases of evolution of massive stars. The evolution 
of stars with radiation-dominated envelopes can require pro- 
hibitively short timesteps in MESA star if the standard mixing 
length theory is adopted. This problem usually appears dur- 
ing the evolution of high mass and/or high metallicity stars af- 
ter hydrogen-core burning and prevents evolution to core col- 
lapse. We discuss in 7.2 our treatment of superadiabatic con- 
vection in these envelopes, which allows uninterrupted evolu- 
tion, from ZAMS to core collapse. 

Since it is relevant to later discussions we start with a plot 
of the OPAL opacity data (Iglesias & Rogers 1996) and 60 Mq 
ZAMS models in Figure 38. The plot is inspired by Figure 1 
of Cantiello et al. (2009). The left-hand panel of Figure 38 
shows the OPAL data for five different Z values at constant 
X = 0.7 and \og(p/T(,^) = -5, where T(, is the temperature in 
units of 10* K. The right-hand panel shows the opacity pro- 
files of five 60 M0 ZAMS models for the same five Z values. 
The model profiles exhibit the same general behavior in the 
opacity-temperature profile as the raw opacity data. Of par- 
ticular importance are the iron opacity bumps that occur at 
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200 kms ' (15MAG). The solid black line shows MESA star results, and the dashed gold line shows the STERN calculations. 
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Figure 29. Same as Fig. 28. As function of mass coordinate we plot the values of the diffusion coefficient for convection (MLT), Eddington-Sweet circulation 
(ES), magnetic torques by dynamo generated fields (TS), Dynamical Shear (DSI), Secular Shear (SSI) and Goldreich-Schubert-Fricke (GSF) instability. Following 
STERN, we turn off the Solberg-H0iland instability (SH) for this comparison. This does not affect the results, as the diffusion coefficient for SH is usually smaller 
than the ones for ES and TS. The values of the specific angular momentum j and the angular velocity Li> are also plotted. Left panel shows the results using 
MESA star, while the right panel shows analogous STERN calculations. 



log T X 5.3 and 6.3. These bumps cause both the local radia- 
tion pressure to dominate and the luminosity to approach the 
Eddington luminosity Z-Edd- 

Where both the pressure is dominated by radiation and Liad 
approaches Ledd, specific conditions can be reached that cause 
convection and inversions in density and gas pressure. To de- 
fine the conditions under which these occur, we follow the 
discussion of Joss et al. (1973), going from high to low L,aci. 
We assume that dT/dr < 0, dP/dr < 0, and that the inertial 
terms in the momentum equation are small. First, we establish 
a condition for the occurrence of an inversion in the gas pres- 
sure i|as- Recasting the equation for the temperature gradient 
gives 



and using the equation of hydrostatic equilibrium, one obtains 

(22) 



d^ad _ ^rad 
Lmd 



dP 



Writing di|a,s/dr = d(P - i?ad)/d'' and using equation (22) and 
the fact that both Z?ad and P monotonically decrease with r, 
one obtains 



dP,. 



dr 



d/?ac 

dr 



^Edd 



(23) 



^rad — 



47Tr^c d7?ad 



PK 



dr 



(21) 



Since diJad/dr < 0, equation (23) implies that for L^ad > ^-Edd^ 
the gas pressure gradient will increase outward, di|as/d?" > 0, 
as shown by Joss et al. (1973). 

The next step is to establish the condition for a density in- 
version to occur. Writing the gas equation of state as /|as - 
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Figure 30. IVIagnetic field structure and angular momentum distribution for model 15MAG at different evolutionary stages (see Table 4). The curves show 
profiles for specific angular momentum (J), angular velocity (ai), azimuthal and radial components of magnetic field (B^ and B,-). The shaded regions represent 
convective parts of the star. Compare with Fig. 1 in Heger et al. (2005). 




Figure 31. Specific angular momentum distribution at dift'erent evolutionary stage for 15MAG and 15ROT. See Table 4 for the definitions of these times. 
Compare with Fig. 2 in Heger et al. (2005). 
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Table 4 

Evolution of Angular Momentum at Fiducial Mass Coordinates for a Z = 0.02, 15 Mq star initially rotating with Ueq = 200 km s"' with and without the inclusion 



of magnetic fields. 
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NOTE: " Results from Table 1 of Heger et al. (2005); * See e.g., Petrovic et al. (2005); Yoon & Langer (2005); Yoon et al. (2006); ^40 % central hydrogen mass 
fraction; ''l % hydrogen left in the core; ""l % helium burnt; ^^50% central helium mass fraction; ^1% helium left in the core; ''central temperature of 5x10* K; 
'central temperature of 1.2x10'' K; ^central oxygen mass fraction drops below 5 %; ^central Si mass fraction drops below 10"'*; 
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Figure 33. Evolution of internal specific angular momentum for the 25MAG model. Solid lines show MESA star result, while dashed lines refer to STERN. 
Left panel shows the evolution from zero age main sequence to He ignition. Right panel shows the evolution during core He-burning (from 50% of He in the core 
to He depletion). Notice the different axis range in the two plots. 
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Figure 34. Evolution to Si-depletion of the internal specific angular mo- 
mentum for the 25MAG model. 



PgAp, Pmd) gives 



dr 



dp 



dp /5Pgas\ di?-ad 



dr 



(24) 



Solving eq. (24) for dp/dr and using eq. (23) eliminates 
dPgas/dr. Gas equations of state have {dPg3s/dp)T > 0, so that 
for dp/dr > (a density inversion), one must have 



dr 



^rad 



- 1 - 



5i?-ad 



> 0. 



Recognizing that d/?ad/dr < 0, we find that a density inversion 



occurs when 



-^rad 



--Edd 



--Edd 



dp. 



rad 



(25) 



This equation is identical to eq. (8) of Joss et al. (1973). Since 
under conditions of interest {dP^a?,ldP^ad)p > 0, we have Li„v < 
^Edd- For Linv < i-rad < L^dd, ^ density inversion will occur 
even though di|as/dr < 0. 

Next, we shall consider the luminosity Lonset at which con- 
vection occurs. In a convective region, the entropy is either 
constant or declining with radius. Hence, convection will oc- 



dln/;-ad ^ / glnj^ad 
dlnP ^ \ dlnP 



(26) 



using equation (22) and solving for the luminosity, we find 
that convection starts once 



^rad ^ ^onset ■ -j^ 

^'Edd i'Edd 



P 

-fga,s 



d In i?ad 
dlnP 



(27) 



Equation (27) corresponds to eq. (9) of Joss et al. (1973). As 
argued in that paper, entropy decreases as density increases; 
therefore a density inversion implies a superadiabatic gradi- 
ent, and as a result, Lonset < ^^inv This can be shown explicitly 
for a chemically homogenous mixture of an ideal gas and ra- 
diation. For such a mixture, equation (25) becomes 



^ad ^ ^inv 
Lsdd L^dd 



l-P.aJP 



1 -3PgaJ4P 



and equation (27) becomes 

^onset 

8(l-i|as//')(4-3Pgas/P) 



^rad 
^'Edd 



32 - 24P„,JP + 3iP.,JP)^ 



(28) 



(29) 



allowing one to show that Lonset < Linv At high lumi- 
nosities where the gas becomes radiation-dominated, how- 
ever, the diff'erence between Lonset and Linv becomes small. 
Expanding equations (28) and (29) for Pg^s/P 1 gives 



-L„ 



(3/4)x(7|as/^')xLEdd- For such high-luminosity. 
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Figure 35. Specific angular momentum distribution for tlie two 16 M© models. In the top panels, the solid curves show the distribution of specific angular 
momentum at dilferent evolutionary stages. The other curves in the top panel show the specific angular momentum of the last stable orbit around a Schwarzschild 
black hole, a maximally rotating Kerr black hole (a = 1 ), and a black hole with a Kerr parameter corresponding to the angular momentum content of the stellar 
progenitor at that mass coordinate. The bottom panels show the free-fall time at the relative mass coordinate at the end of Si-burning. Notice the dift'erent ranges 
of the y-axis. These models can be compared to the calculations of Woosley & Heger (2006), in particular their models 16SG and 16TI respectively. 
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Figure 36. Evolution in the HR diagram for two rotating 40 M© models at 
Z = 10"''. The slower rotating model evolves toward the red part of the HRD; 
the other model evolves toward the blue part of the HRD. The internal evolu- 
tion of the angular momentum is shown in Fig. 37. This can be compared to 
Fig. 2 of Yoon & Langer (2005). 



radiation-dominated stars, a small inefficiency in convection 
is sufficient to drive a density inversion. 

We now demonstrate that such inefficient convection can 
arise in the convective, radiation-dominated, envelopes of 
massive stars. In order of magnitude the convective and ra- 
diative fluxes are, respectively, /Vonv ~ pc, (Vj- - Vaj)"^^^ and 
^rad ~ ci?ad/T. To Carry the flux, we need Fcom ~ ^radi equat- 
ing and substituting pc^ ~ P ~ P„^^, we arrive at an expression 
that sets the level of superadiabaticity. 



(30) 



Under typical conditions in massive star envelopes, c/cs ~ 
10"^ at the iron opacity bump, but at this location, r is not 
large enough to prevent the superadiabaticity from triggering 
a density inversion. 

The lines in Fig. 39 show these luminosity conditions as 
a function of Pgas/P, and reveal that as the stellar conditions 
become radiation dominated, there is only a small gap be- 
tween a convective model that is adiabatically stratified and a 
model with a density inversion. This corresponds to the re- 
gion between the curves Lrad = ^onset (dot-dashed line) and 
Lrad - Lim (dashed line). The gas pressure does not invert un- 
til L > L^dd, which in Fig. 39 is the region above the solid hor- 
izontal line. We show profiles from a 30 Mq (left panel) and 
a 70 Mo model (right panel). These are from the first cross- 
ing of the Hertzsprung gap when T^ff - 5000 K. Each dot 
corresponds to a zone in the calculation; as the profile moves 
outward from center to surface the traces go from bottom to 
top in the plot. The blue dots indicate zones where the star is 
radiative; red indicates convection; a black border denotes a 
density inversion, dp/dr > 0; and yellow indicates a gas pres- 
sure inversion, dPg^^s/dr > 0. There is excellent agreement 
between the detailed MESA evolutionary calculations and the 
analytical conditions (eq. [28] and [29]). The 70 M© profile 
goes into the low P„as/P, high Lrad/i'Edd regime. 

Figure 40 displays the physical conditions in the 70 Mq 
model where the density and gas pressure inversions develop. 
The panels display, from top to bottom, density, gas pres- 
sure, total pressure, and entropy, all as functions of radius. 
The total radius is R = 1330 Rq. Regions with V > Vad and 
Aad < Lim < ^Edd ^rc marked with a small red dot. Regions 
where Liny < i^ad < ^Edd (cf eq. [28]) are marked with a 
large red dot with a black border. Regions where Lrad > ^-Edd 
are marked with a large yellow dot with black border. Al- 
though the pressure (panel c) is well-behaved in this supera- 
diabatic (panel d) region, a density inversion does develop 
where L^d > Aad > Lim (panel a) and a gas pressure in- 
version develops (panel b) where Lrad > L^dd, as predicted. In 
this region the superadiabaticity Vj- - Vad > 10"^ and is greater 
than unity for r/R© > 1300. This is much larger than a typical 
value (~ 10"^) where convection is efficient and results in the 
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are showing the distribution of specific angular momentum at dift'erent stages of the evolution, together with the specific angular momentum of the last stable 
orbit around a Schwarzschild black hole, a maximally rotating Kerr black hole (a = I ) and a black hole with a Kerr parameter corresponding to the angular 
momentum content of the stellar progenitor at that mass coordinate. Note the different ranges of the y-axis. The evolutionary tracks for these models are shown 
in Fig. 36. These calculations should be compared to Fig. 5 of Yoon & Langer (2005). 




entropy decreasing with r as shown in panel c. 

7.2. Treatment of Superadiabatic Convection in 
Radiation-Dominated Regions 

In MESA star the superadiabatic gradient arising in 
radiation-dominated envelopes can force the adoption of pro- 
hibitively short timesteps. Energy is mostly transported by 
radiation, and the convective velocities resulting from MLT 
approach the sound speed. The stability of such radiation- 
dominated envelopes has been questioned in the past, and is 
still a matter of debate (see, e.g., Langer 1997; Bisnovatyi- 
Kogan & Dorodnitsyn 1999; Maeder2009). Hydrodynamical 
instabilities and the transport of energy from waves excited 
by near-sonic turbulent convection are important for energy 
transport, and three-dimensional hydrodynamical calculations 
are required to capture fully the complex physics occurring in 
this regime. 



Here we develop a treatment that allows MESA star to cal- 
culate models of massive stars up to core collapse. For every 
model, MESA star computes the values of 

.Imax = max^^^j and ySmin = min . (31) 

When jSjnin is small and A^^,^ is large, and MLT yields a <5v > 
^v.thresh, we artificially decrease the superadiabaticity, 5y = 
- Vad, implied by MLT. The default of the user-specified 
parameter 5v,thresh is sufficiently large, ~ lO"'', much larger 
than the typical values for efficient convection. 

MESA star sets Vj- to reduce the 5v ~ ^v,thre,sh by ^ factor 
Q'v/v, where fy is specified by the user, and ccv is updated at 
each timestep to a linear combination of its previous value and 
a value av('imax,y6min)- For large values of /Imax and small val- 
ues of jStnin, ov — » 1 ; in typical usage, the transition happens 
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Figure 39. The critical luminosities L,-ad = ^onsct (eq- [29], dot-dashed line), 
A-ad = Anv [28], dashed line), and Li^j = Ledd (solid line) as a func- 
tion of Pgas/P for an ideal gas-radiation mixture. Compare this with Fig. 1 
of Joss et al. (1973). For Lraj < Lon.sct> the gas is convectively stable; for 
^onset < irad < ^-inv the gas is convective; for Li„v < ^rad < L^dd, the den- 
sity is inverted, dp/dr > 0; and for L^dd < Aad> the gas pressure is inverted, 
d^as/dr > 0. Overlaid on the plots are the profiles from a 30 Mq (left panel) 
and a 70 M© (right panel) model with Z = 0.O2: blue dots indicate zones that 
are radiative; red dots indicate Vjaj > V;,j; dots with a black border have a 
density inversion; and the yellow dots with black borders indicate a gas pres- 
sure inversion. As the profile moves out from the stellar center it traces out 
the points on the plot from bottom to top. Only a part of the model profiles 
are visible in the plot. The calculations correspond to the first crossing of the 
Hertzsprung gap when T^g = 5000 K. 



where A^^^ ^ 0.5 and /3jnm ~ 0.3. For small values of A^^^ 
and large values of jSmin, 5v ~> 0. Thus fy sets the maximum 
reduction of (Jv - ^v.thi-esh- Figure 41 shows how IlESA star 
turns on the reduction in 5v as a star evolves. Tracks in the 
HR diagram are shown for four stellar models: 15, 25, 30, 
and 70 Mq. The color of each line indicates the value of av at 
each point. 

Such a decrease of the temperature gradient reduces Ljad 
and implies additional physical transport. Potential agents 
for the excess transport include waves excited by turbulent 
convection (see, e.g., Maeder 1987a) and radiative diffusion 
enhanced by porous clumping of the envelope (e.g., Owocki 
et al. 2004). As these radiation-dominated envelopes might 
be physically unstable, with a resulting strong enhancement 
of mass loss, we caution that the results of any ID stellar evo- 
lution calculation for the late evolutionary phases of massive 
stars should be considered highly uncertain. 

We now show a comparison of MESA star calculations of 
rotating massive stars done with and without MLT-n-. We 
used the 25 Mg model described in §6.4, which at Z = 0.02 
is around the upper mass limit that can converge using a rea- 
sonably short timestep without having to rely on the IVILT-n-. 
The most prominent difference between the calculations is the 
evolutionary track in the HR-diagram (Fig. 42). This is not 
surprising, as IVILT and JVILT-fH- result in different efficiencies 
of energy transport in radiation-dominated stellar envelopes. 
The sharp drop in L for the MLT-I-+ case is the result of a 
brief period of enhanced mass loss due to super-critical rota- 
tion. The structure and the angular momentum content of the 
collapsing core are weakly dependent, however, on the choice 
of IVlLT vs. MLT++ (Fig. 43). 

7.3. Core-Collapse Progenitor Models 




Figure 40. The panels display, from top to bottom, the density, gas pres- 
sure, total pressure, and entropy as functions of radius for the 70 M© model 
shown in Fig. 39. The range of radii is restricted to the region where density 
and gas pressure inversions develop. Each zone is marked by a dot; a small 
red dot indicates convection with no predicted gas or gas pressure inversion 
(Aad < Anv); a large red dot with black border indicates a predicted density 
inversion but no gas pressure inversion (Li„v < Lmd < ^Edd); and a yellow dot 
with black border indicates a convective region with a predicted gas pressure 
inversion (Liad > ^Edd). The total pressure (panel c) is well-behaved at all 
radii. Note also the decrease in entropy (panel d): the region is superadia- 
batic. 
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Figure 41. HR diagram of 15, 25, 30, and 70 Mq models. The color indicates 
the value of av at that point in the star's evolution. For the 25 Mq and 30 Mq 
stars, there is a sharp spike in (j-y as the star crosses the Hertzsprung gap 
followed by a sharp drop at the base of the red giant branch. The 70 Mq 
model has ay > 0.9 for its entire evolution. 
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black lines show MESA star results with MLT (black) and MLT++ (orange), while the dashed blue line refers to STERN calculations. The star symbol shows 
the location where we started the calculation for the RSG pulsations discussed in §7.4. 

et al. (2009). For Teff > lO'^K and H-surface fraction > 0.4, 
the mass-loss prescription of Vink et al. (2001) is used. In 
the same temperature range, but when the H-surface fraction 
decreases below 0.4, Nugis & Lamers (2000) determine the 
mass-loss rate. At low temperatures (Tgff < 10"^ K) the mass- 
loss rate of de Jager et al. (1988) is used. 

Figure 44 shows the central conditions of these massive ro- 
tating models. For each model the calculation stops when 
any part of the collapsing core reaches an infall velocity of 
1000 km s"' . Some of the initial and final properties are sum- 
marized in Table 5. These calculations are performed to reveal 
the new capabilities of MESA star. The values of the param- 
eters for these calculations have not been calibrated against 
existing calculations or observations. 
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Figure 43. Evolution of internal specific angular momentum for the two 
25 Mq models of Fig. 42. The dashed lines show models calculated with 
MLT++. Due to dift'erent excursions in the HRD (see Fig. 28) calculations 
with MLT and MLT++ end with different final masses. There are no sub- 
stantial changes, however, in the specific angular momentum content of the 
stellai' cores. 

We evolve a grid of massive stars initially rotating with 
O/Qciit = 0.2. The models have been initialized using solid 
body rotation. Models with initial M/Mq = 30, 40, 50, 60, 
70, 80, 90 and 100 have initial Z = 0.02, while models with 
initial M/Mq = 120, 150, 250, 500 and 1000 have been ini- 
tialized with Z = 0.001. To calculate convective boundaries 
we adopt the Ledoux criterion including the impact of semi- 
convection (with Use - 0.02, see §4.1). The transport of an- 
gular momentum and chemicals by rotational instabilities and 
magnetic torques is included and calibrated following Heger 
et al. (2000, 2005) and Yoon & Langer (2005). Wind mass- 
loss is been implemented following the recipe of Glebbeek 



7.4. Radial Instability of Red Supergiants 

Massive red supergiants (RSG) are unstable to radial pul- 
sations driven by the A--mechanism in the hydrogen ionization 
zone. Both linear and non-linear calculations show the oc- 
currence of oscillations with the period and growth rate of 
the dominant fundamental mode increasing with L/M (Li & 
Gong 1994; Heger et al. 1997; Yoon & Cantiello 2010). The 
periods are of the order of years. As discussed by Yoon & 
Cantiello (2010) the occurrence of RSG pulsations can im- 
pact stellar mass-loss rates and modify the evolution of mas- 
sive stars above a certain mass. We study the occurrence of 
RSG pulsations with MESA star and compare results with ex- 
isting non-linear calculations. 

In Fig. 45 we show the capability of MESA star to exhibit 
radial oscillations in luminous RSGs. We use the same 25 Mq 
rotating model discussed in § 6.4, and we restart the calcu- 
lation when the He mass fraction in the core is Yc = 0.7. 
For non-rotating RSG with Z - 0.02, Yoon & Cantiello 
(2010) found pulsation periods in the range 1-8 yr. To re- 
solve the RSG pulsations we force the timestep to < 0.01 yr, 
much shorter than the usual timestep during He-burning {6t > 
10^ yr, see Appendix B.3). This explains why RSG pulsa- 
tions are usually not found during the evolution of massive 
stars. Before the code stops due to the emergence of super- 
sonic radial velocities in the envelope, we find a pulsational 
period * 4 yr, in good agreement with the results of Yoon & 
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Table 5 

Initial (ZAMS) and final (core-collapse) properties. 





Ziiii 


i2 / ficrit 








Arn^ 


A^He" 


m/ 






7Fe' 


[Mo] 






[kms-'] 


[ergs s] 


[Myr] 


[Myr] 


[Myr] 


[Mq] 


[Mq] 


[ergs s] 


[ergs s] 


30 


0.020 


0.20 


129.69 


3.28x10" 


6.30 


5.87 


0.36 


17.77 


1.41 


2.87x1050 


1.03x10"^ 


40 


0.020 


0.20 


122.86 


4.87x10" 


5.06 


4.71 


0.31 


19.37 


1.81 


3.77x10^0 


LeixlO""* 


50 


0.020 


0.20 


112.02 


6.30x10" 


4.41 


4.08 


0.29 


25.04 


1.38 


5.39x10^0 


1.09x10""* 


60 


0.020 


0.20 


98.37 


7.34x10" 


4.04 


3.66 


0.35 


22.88 


1.76 


7.81x10^° 


2.76x10""* 
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78.76 
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26.19 


1.75 
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1.54x10""* 
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3.38 


0.29 
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3.10 
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0.27 
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0.26 
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4.79x1 0^' 
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2.09x10" 


1.99 


1.77 
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" Initial rotation rate, see definition in §. 6.3. 
''Initial equatorial rotational velocity. 
■^Total initial angular momentum, 
''stellar lifetime. 

■^Main sequence and core He-burning lifetimes. These are defined as the interval between onset of core burning and depletion of central hydrogen (or helium) 
to 1 % by mass. 

f- 



Final mass. 
^Mass of the Iron core (if present). 

Final total angular momentum. 
'Final total angular momentum of the iron-core, 
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Figure 44. Evolution of and pc in the massive rotating models. The 
locations of core helium, carbon, neon, oxygen, and silicon burning are la- 
beled. A dashed curve marks the electron-positron pair-instability region 
where Fi < 4/3. All models are rotating initially at 20% of critical rotation. 
The calculations include the eft'ects of rotation and Spruit-Tayler magnetic 
fields as discussed in §6. Models with initial mass < 100 M© have initial 
metallicity Z = 0.02, while models with mass > I2OM0 have initial metal- 
licity Z = 0.001. The end of the line for each mass corresponds to the time 
of core-collapse, defined as when any part of the collapsing-core exceeds an 
in-fall velocity of 1 000 km s" ' . The tracks for the 80 Mq and 90 Mq overlap 
in this plot. 

Cantiello (2010). 

8. SUMMARY AND CONCLUSIONS 
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Figure 45. Surface properties of a pulsating RSG. This is the same 25Mo 
model discussed in Sec. 6.4, evolved from fo = 6.851 Myr (corresponding 
to F(. = 0.7, star symbol in Fig. 32) with timesteps St < 0.01 yr. The black 
line shows the evolution of the stellar radius, while the orange line shows the 
value of the surface radial velocity (in units of the local sound speed). The 
inset shows the corresponding evolution in the HR-diagram. 



We have explained and, where possible, verified the im- 
provements and major new capabilities implemented in MESA 
since the publication of Paper 1. These advancements include 
evolutionary modeling for giant planets (§2), tools for astero- 
seismology (§3), implementation of composition gradients in 
stellar interiors and their impact on convective mixing (§4), 
the evolution of intermediate mass stars and white dwarfs 



Modules for Experiments in Stellar Astrophysics (MESA) 



29 



(§5) the treatment of rotation during stellar evolution (§6), 
addressing the onset of radiation pressure dominance in the 
envelopes of evolving massive stars due to the iron opacity 
bump, and evolving massive stars to the onset of core collapse 
(§7). The enhancements include the physics modules (Ap- 
pendix A), the algorithms (Appendix B), and the addition of 
a MESA Software Development Kit (Appendix C). MESA star 
input files and related materials for all the figures are avaliable 
at http : //raesastar . org. 

These hitherto unpublished advancements have already en- 
abled a number of studies in planets (e.g., Passy et al. 2012; 
Huang & Gumming 2012; Carlberg et al. 2012), classical no- 
vae (Denissenkov et al. 2013), asteroseismology (e.g., Yang 
et al. 2012; Burkart et al. 2012b; Moravveji et al. 2012), rota- 
tionally induced mixing (e.g., Denissenkov 2010; Chatzopou- 
los et al. 2012; Chatzopoulos & Wheeler 2012) and enabled 
the discovery of new features in the evolution of low-mass 
stars (Denissenkov 2012). In addition, these enhanced ca- 
pabilities have allowed for applications of MESA star that 
were not initially envisioned, such as explorations of stars un- 
der modified gravity (Chang & Hui 2011; Davis et al. 2012), 
and stellar oscillations induced by tidal disturbances in dou- 
ble white dwarf binaries (Fuller & Lai 2012a,b; Burkart et al. 
2012a). 

As an open source "instrument" for stellar astrophysics, it 
is difficult to predict all the ways in which future develop- 
ment of MESA will occur. We do know, however, that future 
versions of MESA will include advances in physics modules, 
features driven by the MESA user community, and architec- 
tural refinements. For example, the plethora of asteroseis- 
mological data is driving new initiatives to incorporate non- 
adiabatic pulsation codes, where possible, into MESA. The 
prevalence of interacting binary star systems, especially for 
massive stars, has increased the pressure for MESA develop- 
ment efforts that would yield the capability to simultaneously 
evolve two interacting stellar models. Physics module devel- 
opments will likely include general relativistic corrections to 
the stellar structure equations (e.g., difference between gravi- 
tational and baryonic mass), the mass diffusion coefficients in 
electron degenerate environments, phase separation in cooling 
white dwarfs, and nuclear statistical equilibrium solvers. We 
also expect the transition from multicore systems (with order 
10 cores) to many-core architectures (with order 100 cores) to 
drive new directions in MESA's algorithmic and architectural 
development. 
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APPENDIX 

UPDATES TO INPUT PHYSICS MODULES 

There have been many updates and improvements to the physics modules since Paper I. In this appendix, we describe the 
changes that have been made to the microphysics modules chem (§ A.l), eos (§ A. 2), kap (§ A. 3), and net (§ A.4). We conclude 
by listing updates to the atmosphere boundary conditions (§ A. 5). 

Atomic and Nuclear Data 

The chem module now has the latest version (v2.0) of the JINA reaclib nuclide data (Cyburt et al. 2010). This contains 
updated mass evaluations, and now includes 7853 nuclides up to ^^^Cn. For precision work, the chem module now distinguishes 
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between the atomic mass number A, — the number of nucleons in a given isotope — and the atomic mass W,-. The abundance of a 
species / is defined as 

ri: 

Y, = —, (Al) 
where «b is the baryonic number density. The baryon fraction Xj is then 

X, = y,A, = — , (A2) 

"B 

Note that 2/^/ - "b/wb = 1 and is invariant under nuclear reactions. We then define the baryon density (in mass units) as 

P = wb^u, (A3) 

where /«□ = 1.660538782 x lO"^'* g is the atomic mass unit (CODATA 2006 value; Mohr et al. 2008). Note that the numerical 
value nia, along with other physics constants, are defined in the const module. The atomic mass of isotope / is defined in MESA 

as 

W, = A,- + ^ , (A4) 

where A,7c^ is the mass excess of isotope /. This treatment neglects the electronic binding energy, and A is therefore independent 
of the ionization state of a given species. The electron rest masses are, however, included in this definition, since the W, are 
atomic masses. 

The MESA micrqphysics modules — kap, eos, neu, and net — use p, T, and {Xj} as inputs. MESA star multiplies p by a mass 
correction factor W/A = 2, W,y,7 2, A,y, to distinguishes between A, and W, before starting the calculation for a timestep. A 
call to the routine composition_info in the chem module returns the following averaged quantities: the mean atomic mass 
number, A = 2, y,A,/ 2, F,, mean atomic charge number, Z = 2, Z,y,7 2,- F,-, mean square atomic charge number, 2, ZfYj/ 2, Yi, 
the electron abundance, = Z/A, and the mass correction term, W/A. In addition, the routine returns the derivatives of A, Z, 
and W/A with respect to the baryon fractions X,; 



dA 

dXi 

5Z 
dXi 

d{W/A) 



dXi 



1 



Wi_ _w\ 

Ai ~ Jj l.iXi 



(A7) 



Note that the routine does not make any assumption in these derivatives that 2, X, = Y^iAiYi - 1; in this formulation, 2,- X, is not 
explicitly set to unity. 

At the beginning of each Newton iteration, the abundances are checked. A mass fraction is considered good if its value exceeds 
min_xa_hard_liinit. If all mass fractions meet this standard, then the mass fractions are clipped to range from to 1, and the 
mass fractions are summed. If the sum differs from unity by less than a value sum_xa_tolerance, then the mass fractions are 
renormalized to sum to unity; otherwise, the code reports an error. Currently composition derivatives are ignored in the eos and 
kap routines. Equations (A5)-(A7) allow, however, future additions to these routines to compute these derivatives analytically. 

Equation of State 

The only significant change to the eos module since Paper I is the addition of tables for Z > 0.04, where Z is the mass fraction 
of all elements heavier than He. The eos module as described in Paper I supplied equation of state (EOS) tables for Z - 0.0, 
0.02, and 0.04 at temperatures and densities for which neutral and partially-ionized species are present (see Paper I, Figure 1). 
For Z > 0.04 MESA switched to the HELM EOS (Timmes & Swesty 2000), which assumes full ionization. In order to rectify the 
inconsistent treatment of the partially-ionized region at high Z, new EOS tables have been computed (J. MacDonald, priv. comm.) 
using the MacDonald EOS code (MacDonald & Mullan 2012) for Z = 0.2 (scaled-solar), and two Z = 1.0 compositions: one 
with 49.5% C, 49.5% O, and 1% scaled-solar by mass; and one with 50% C and 50% O by mass. Here "scaled-solar" refers to 
the Grevesse & Noels (1993) solar heavy element distribution adopted in the OPAL EOS tables (Rogers & Nayfonov 2002). 

Opacities 

The kap module now divides the opacity tables into a high-temperature domain, log(r/K) > 4, and a low-temperature domain, 
log(r/K) < 4; the exact range of log T over which the tables are blended can be adjusted at runtime. This treatment differs 
from the opacity tables described in Paper I, which combined high- and low-temperature opacities into a single set of tables. The 
motivation for separating the tables is to facilitate using different sources of low-T opacity data. The kap module now supports 
low-T opacities from either Ferguson et al. (2005) or Freedman et al. (2008) with updates to the molecular hydrogen pressure- 
induced opacity (Frommhold et al. 2010) and the ammonia opacity (Yurchenko et al. 2011). Either set may be selected at run 
time. The electron conduction opacity tables, based on Cassisi et al. (2007), have been expanded (Potekhin 2011, priv. comm.) 
to cover higher temperatures (up to 10'" K, originally 10^ K) and densities (up to 10"^ g cm"^, originally 10"^^^ g cm"^). 
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Nuclear Reactions 

Substantial improvements to the net module have been made since Paper 1 to increase the flexibility of the nuclear reaction 
networks (see § B for working details). One such improvement is the standalone one-zone burn routines. These now operate on 
a user-defined initial composition, nuclear network, and a thermodynamic trajectory. Choices for the thermodynamic trajectory 
include a burn with density and temperature held fixed, a burn with pressure held fixed, and a burn with the density and temper- 
ature following an arbitrary, user-specified profile. This last option is activated by setting read_T_RhoJiistory= . true . and 
specifying the file name containing the profile through the variable TJlho Jiistory_f ilename. The MESA one-zone burn routines 
now include user-specified options for the family of stiff ordinary differential equation integrators from Hairer & Wanner (1996). 
In addition, three user-defined switches are provided to switch between using dense matrix linear algebra solvers, for smaller 
networks, and sparse matrix linear algebra solvers for larger ones. The option decsol_switch sets the number of isotopes at 
which the switch occurs; options small_mtx_decsol and large_mtx_decsol specify the dense and sparse solvers, respectively. 

Figure 46 shows the constant pressure option of these routines operating on conditions that might be encountered for helium 
burning on the surface of a white dwarf. The initial pressure is 3.1 x 10^^ ergs cm"', the initial temperature is 2 x 10^ K, the initial 
composition is Xi'^Ht) - 0.98 and X('^N) = 0.02, and the system was evolved for 10^ s with a 19-isotope network. Evolution 
of the density and temperature under the constant burn conditions are shown in the lower panel of Figure 46. The temperature 
slowly increases and the density slowly decreases as the material begins to burn and release energy at a rate of tnuc - cp dT jdt. 
When the temperature crosses a critical threshold at a; 20 s, a runaway occurs as the temperature rapidly rises and the composition 
bums to heavier elements. The material then establishes a final equilibrium state, no energy from nuclear burning is injected into 
the system, and the temperature reaches a plateau. 

The upper panel of Figure 46 compares the evolution of key isotopes and the energy generation rate per unit mass of the 
MESA one-zone burner (colored and labeled curves) with an independent one-zone burner (dashed black curves) based on Timmes 
(1999). These comparisons indicate that both one-zone burns produce a final composition that is mostly ^''Ti and "^^Cr. Over most 
of the evolution, the two calculations give mass fractions of various isotopes that agree to within 2-3 significant digits. Larger 
differences in some of the heavier isotopes at the end of the calculation are due to differences in the adopted nuclear reaction 
rates. 




Time (s) 

Figure 46. A one-zone helium and nitrogen bum at constant pressure, P = 3.1 X 10^- ergs cm"^, starting from an initial temperature of T = 2x 10** K. Evolution 
of the temperature and density are shown in the lower panel, while the upper panel shows the mass fraction of key isotopes (right axis) and the energy generation 
rate per unit mass (left axis; red curve). MESA results are shown by the colored and labelled curves, and the results from an independent one-zone burner (Timmes 
1999) are shown by the dashed black curves. 



Another improvement is the net module now accesses reactions from both weaklib and reaclib. Rather than evaluating 
the standard seven-parameter fit for Na{o-v} for the reaclib rates (Cyburt et al. 2010) every time a reaction rate is needed, the 
net module caches separate rate tables for each reaction. Inverse rates are calculated directly from the forward rates (those 
with positive Q-value) using detailed balance, rather than using fitted rates. This is important for explosive nucleosynthesis 
approaching nuclear statistical equilibrium (see Calder et al. 2007). The nuclear partition functions used to calculate the inverse 
rates are taken from Rauscher & Thielemann (2000). 
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Atmosphere Boundary Conditions 

The atm module provides the surface boundary condition for the interior model. A collection of four new options that extend 
the set described in Paper I are now available. 

1. solar JIopf_grey: Implements the solar-calibrated Hopf T(t) relation, where 

T\T)^^T^g[T + qiT)], (A8) 

and 

q(T) ^qi+qi exp(-q'3T) + q^ expi-q^T). (A9) 

The qi are fit to the solar atmosphere with resulting values qi = 1.0361, q2 - -0.3134, q^, - 2.448, q^ - -0.2959, and 
qs = 30.0 (J. Christensen-Dalsgaard, 2011, priv. comm.). 

2. grey_and_kap: Expands on the simple assumption that P ^ rg/K by iterating to find a consistent solution among P, T, and 
K(p,T). 

3. grey .irradiated: Implements the Guillot (2010) T(t) relation that includes both external irradiation by the star and 
cooling flux from the interior; see Guillot (2010, eq. 49) along with the discussion and results in §2.4. In addition to the 
external and internal fluxes, this boundary condition requires two constant opacity values: for the external radiation, and 

for the thermal radiation generated within the atmosphere. This boundary condition is unique in that it is applied at a 
specified pressure level, as opposed to optical depth. This pressure must be chosen sufficiently high to capture any heating 
of the atmosphere by the irradiation. 

4. WD_tau_25_tables: Provides as outer boundary conditions the values of 7|as and T at log(T) = 1.4 as extracted from 
pure hydrogen model atmospheres of WDs (Rohrmann et al. 2012; Rohrmann 2001). The tables span a range of effective 
temperatures and surface gravities: 2, 000 K < T^ff < 40, 000 K and 5.5 < \ogg < 9.5. See § 5.2 for an example of the use 
of these tables. 

NUTS AND BOLTS 

We now briefly describe the primary components of evolution calculations. MESA star first reads the input files and initializes 
the physics modules to create a nuclear reaction network and access the EOS and opacity data. The specified starting model is 
then loaded into memory and the evolution loop is entered. 

Evolve a Step 

The top level routine for evolving a star for a single timestep is do_evolve_step. If this is the first attempt to do a step 
starting from the current state, the model is remeshed (see §B.4), and information for MLT-H-i- is prepared by the routine 
set_gradT_excess_alpha (see §7.2). Sufficient information is saved so that if necessary it will be possible to make other 
attempts (i.e., after a redo, a retry, or a backup). In addition to the current state, we keep the previous state (called "old"), and the 
one that came before "old" (called "older"). During the step, the current state is modified, and the old one holds the state at the 
start of the step. If we do a redo or a retry, we copy old to current to restore the starting state. If we do a backup, we copy older 
to old before copying old to current, making us start at the state prior to the current one. Note that the duration of the timestep is 
determined before the call on do_evolve_step by the process described in §B.3. 

After remeshing and the other initial preparations, do_evolve_step begins the operations that are done on every attempt. It 
first calls the routine do_winds which sets M based on the current radius, luminosity, mass, metallicity, and other properties as 
needed. During the evalution of do_winds there is a call on the user-defined other_wind routine giving users an easy way to 
define different schemes for setting M. 

Information for evaluating the Lagrangian time derivatives is stored by a call to save_£or_d_dt. The ensuing call to 
do_adjustjTiass adds or removes mass without changing the number of grid points (see §B.5). Information for evaluating the 
Lagrangian time derivatives is updated at this point. Variables for the model are evaluated to reflect the changes made by remesh- 
ing and changing mass. This includes evalution of the Brunt-Vaisala frequency (see §3.3), and the diffusion coefficients for the 
mixing of composition (see §4.1 and 4.2). The user-defined routines other_brunt, other jmlt, and other jnixing are called as 
part of this. If rotation is enabled, there is a call to set_rotation_mixing_info (see §6) which in turn calls other_am_mixing. 
If element diffusion from gravitational settling and chemical diffusion is active, the routine do_element_dif£usion adjusts the 
composition and includes a call on other_diffusion The ensuing call to do_struct_burn_mix solves for the new structure 
and composition of the star through repeated Newton iterations (see §B.2). Non-convergence causes do_evolve_step to return 
with a result indicating a failure. Convergence is followed by a call to the routine do_solve_omega_mix which adjusts the total 
angular momentum by solving a diffusion equation (see §B.6); it calls other_torque. There is an option to repeat the operations 
described in this paragraph in case rotationally enhanced mass loss (see §6.3) has not been sufficient to eliminate super-critical 
surface velocities. In such a situation, the mass loss is adjusted iteratively until slightly sub-critical velocities result. In effect, 
this is an implicit solution for the appropriate M when super-critical rotation occurs. 

Next, if specified by the user, smooth_convective_bdy is called to smooth abundances behind retreating convection bound- 
aries. Finally, a call to do_report gathers information and metrics about the timestep for the user. This information will then be 
available to the user's extras_checkjnodel routine. 
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Solving the Coupled Structure, Burn, and Mix Equations 

A call to do_struct_burnjnix invokes a Newton method — an N-dimensional root find — to solve a system of nonlinear 
differential-algebraic equations for the new structure and composition of the stellar model. Here is the number of zones in the 
current model times the number of basic variables per zone and can exceed 100,000. 

The equations to be solved are written as the relation F(basic_vars) - 0, where F is the vector-valued function of the 
residuals. If the basic_vars were a perfect solution to the equations, we would have F = 0; in practice, the solution is never 
perfect. The solution strategy is to iteratively adjust the values of the basic_vars to reduce F towards zero. An approximate 
solution is accepted depending on both the magnitude of F and the relative size of the adjustments to basic_vars. Adjustments 
are chosen using the Jacobian matrix of partial derivatives of all the F equations with respect to all the basic_vars. 

Figure 47 figure shows the three blocks making up the row of the block tridiagonal Jacobian matrix for the tenth from the 
center cell of a 2.5 M© ZAMS model, with black dots showing non-zero entries. The partial derivatives of the equations for cell 
k form the rows of the blocks. In this case, we have 4 equations for the structure of the model (f , T, L, and r) and 8 equations for 
the chemical abundances ('H through ^"^Mg). Each block of the tridiagonal matrix is demarcated by dashed black vertical lines. 
The block matrix on the left shows the dependencies of the equations for cell k on the variables of cell k-1, the one in the middle 
shows the dependencies of the equations for cell k on the variables of cell k, and the one on the right shows the dependencies of 
the equations for cell k on the variables of cell ^ -i- 1 . The dashed lines partition each each block into four sub-blocks to highlight 
the structure and abundance portions of each block. 

The structure of the lower-right subblocks in the left and right blocks shows that the chemical abundance of a particular species 
in cell k depends on the chemical abundances of that species in cells k-l and k+l; this is because of mixing between neighboring 
cells. The lower-right subblock of the center block also shows the interdependencies of abundances due to nuclear reactions in 
the cell. The bottom-left subblocks are zero in the left and right blocks but show dependencies on the P and T variables of the 
center block. This is because the nuclear reactions that change the abundances depend of P and T of that cell but do not depend 
on P and T in the neighboring cells. The columns for L and r are zero in the center lower-left subblock because the equations 
for the abundances do not directly depend on those variables. The upper-right subblocks are zero in the left and right blocks but 
show that the equation for L depends on the abundance variables in the center block. This is because the L equation includes 
results from nuclear burning, and that depends directly on the composition of cell k but not on the composition in neighboring 
cells. The other rows in the center upper-right subblock are zero because the equations for P, T, and r do not directly depend on 
composition. 

Finally, consider the upper-left subblocks that show the dependencies of structure equations on structure variables. The upper- 
left subblock in the center shows that each structure equation in k depends on 3 or 4 of the structure variables in k. The P and 
T equations for cell k also depend on both of the variables P and T in k - I. while the L and r equations for k depend on the 
corresponding variables in ^ - 1 . This pattern reflects the form of the finite differences in the implementation of the structure 
equations: P and T differences use the outer neighbor {k - 1) while L and r differences use the inner neighbor (k + 1). The L 
and r equations for innermost cell k - n use Lcenter and ^center; the P and T equations for the outermost cell k = 1 use the surface 
boundary conditions. 

Variables 



o 



w 
c 
o 

CO 

cr 
LU 







cell k-1 




cell k 




cell k+1 




d{struct|^) 




d(struct|^) 


d(struct|^) 


d{struct|^) 






d{struct|^.-| ) 




d(struct|^) 


d{chem|^) 


d{structk+l) 




p 


• • 




• • 








T 


• • 




• • • •; 








L 










• i 




R 






• • •; 








1H 




• 




»••••• 


:( 


i 


3He 




• 








• 


4He 




• 




»•••••• 




• 


12c 




• 




» • • • 




• 


14n 




• 




» • • • • 




• 


160 




• 




» • • • • 




• 


20Ne 




• 




• • • • 




• 


24Mg 




• 




• • • 




• 






d(chem|^) 


d(chem|^) 


d(chem|^) 




d(chem|^) 






d(chem|^.-| ) 


d(struct|^) 


d(chem|^) 




d(chemk+l) 



Figure 47. One row of the block tridiagonal Jacobian matrix for a 2.5 Mq ZAMS model, with black dots showing the locations of non-zero entries. 



The structure variables for each zone always include the zone average of the natural logarithm of the temperature. In T, the 
luminosity at the outer edge of the zone, L, the natural logarithm of the radius at the outer edge of the zone. In r, and a second 
thermodynamic variable — either the zone average of the natural logarithm of the mass density, Inp, or the the zone average 
of the natural logarithm of the gas pressure, ln/|as- Ideally it would not matter whether Inp or In/^as was used as the second 
thermodynamic variable — for a given temperature and composition the equation of state permits going back and forth between the 
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two. Microphysics packages tend to use mass density as a primary input (i.e., they use a Helmholtz free energy basis) leading to 
the common choice of In p. However, the structure equations are solved only to within a finite but non-zero residual (see above). 
Approximately correct values for the density and temperature can then lead to anomalous pressure profiles, with tiny violations 
of hydrostatic balance. These local violations tend to appear near large jumps in density, such as at a sharp H/He boundary. 
Using i^as as the second thermodynamic variable (effectively using a Gibbs free energy basis) removes these anomalous pressure 
profiles. For example, in stellar models without overshooting or semiconvection, the H/He boundary is extremely sharp. Using 
the gas pressure as the second thermodynamic variable results in single zone step function transitions in the abundances and in the 
density, while the temperature and pressure are smooth across the transition. Applications that demand smooth pressure profiles, 
such as pulsation analysis (see §3), should generally specify the gas pressure as the second thermodynamic variable. 

If the optional hydrodynamic mode is activated, then the radial velocity at the outer edge of the cell, v, is added to the structure 
variables. Figure 47 shows the order of the model variables in the Jacobian: each cell includes the structure variables followed 
by the mass fraction X, of each isotope. Mass and the local angular velocity u are not treated as structure variables because they 
are held constant during the Newton iterations. The mass is set before the iterations, while w is determined after convergence. 
This is computed taking into account loss/gain of angular momentum during the time step, the new stellar structure and internal 
transport of angular momentum calculated by a diff'usion equation (see § B.6). 

The program flow to solve the coupled structure, burning, and mixing equations is to first create the matrix of partial deriva- 
tives using the current candidate solution, solve the block tridiagonal system of linear equations for the corrections to the basic 
variables, apply the possibly damped corrections (see next paragraph) to update the candidate solution, and calculate the residual 
F. Then, if the residual is small enough, we declare victory, otherwise we repeat the general flow with the updated candidate 
solution. 

Each iteration of the Newton solver uses a linear approximation to create a vector of corrections to the model. These corrections 
do not include the physical requirement that the abundance mass fractions need to remain positive. To reduce the possible 
occurrence of negative abundances MESA star now uses a damped Newton scheme. This checks for proposed corrections that 
would produce negative abundances and multiplies the entire correction vector by a factor less than one, so that only part of the 
the full correction is applied. In many cases, this is sufficient to significantly improve the convergence properties of a model. In 
other cases, the damped correction scheme may force so many small corrections that the Newton solver cannot converge within 
the user-specified maximum number of iterations, forcing the previous model to be attempted again with a smaller timestep 
(termed "a backup"). On balance, this is usually a small price to pay for an improved conservation of species and more accurate 
solutions. 

The modules in star provide routines to evaluate the residual equations and create the Jacobian matrix. Given a candidate 
solution (i.e., the set of basic variables for each cell), the microphysics for each cell (EOS, thermal neutrino loss, opacity, nuclear 
reaction rates) are evaluated in parallel (see §B.8). The Jacobian matrix is then further populated with elements from rotation, 
artificial viscosity, and mixing length theory for the temperature gradient, and these are also evaluated in parallel. Each of the 
routines that evaluate these components returns output values and partial derivatives of the output values with respect to the input 
values. Analytic partial derivaties are used whenever feasible, otherwise numerical partial derivatives are supplied. 

Timestep Controls 

Control of the timestep is a critical part of stellar evolution and requires careful trade-offs. The timesteps must be small enough 
to allow convergence in comparatively few iterations but large enough to allow sufficiently efficient evolutions. Changes to the 
timestep must respond rapidly to varying structure or composition conditions, but they need to be controlled to avoid large jumps 
that can reduce the convergence rate or the accuracy of the results. The routine pick_next_timestep performs timestep control 
as a two-stage process. The first stage proposes a new timestep using the H211B low-pass filter (Soderlind & Wang 2006), a 
scheme based on digital control theory. The second stage implements a wide range of tests that can reduce the proposed timestep 
if certain selected properties of the model are changing too much in a single timestep. 

For the first stage, routine hydro_timestep sets the variable for the next timestep, dt_next, according to the relative mag- 
nitude of changes to the basic_vars. The variable reflecting the size of these changes is called varcontrol and is calculated 
by the routine eval_varcontrol. For improved stability and response, the low-pass controller uses previous and current values 
of varcontrol to make the next timestep match the varcontrol_target, v,, which is 10 by default. To make this explicit, 
let (5f,, and (5f,+i be the previous, current, and next timestep, respectively, while Vc,,-i and Vc,, are the previous and current 
values of varcontrol. The maximum timestep for model / + 1 is then determined by 
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where f{x) = 1 H- 2 tan"' [0.5(x - 1)]. This control scheme allows rapid changes in the timestep without undesirable fluctuations. 

The timestep proposed by this low-pass filtering scheme can be reduced in the second stage according to a variety of special 
tests that have hard and soft limits. If a change exceeds its specified hard limit, the current trial solution for the new step is 
rejected, and the code is forced to do a retry or a backup. If a change exceeds its specified soft limit, the next timestep is reduced 
proportionally. The current classes of special cases that can reduce the next timestep are limits based on: 

1 . Number of Newton iterations required to converge. 

2. Maximum absolute change in the mass fraction of hydrogen or helium in any cell. 

3. Maximum relative change in any mass fraction at any cell. 
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4. Magnitude in the relative change in the structure variables in each cell. 

5. Nuclear energy generated in each cell for several categories of nuclear reactions. 

6. Changes in the luminosity resulting from nuclear burning. 

7. Changes at the photosphere in In L and T^s. 

8. Changes in Inpcenter, In Tcenter, ^(H)centei-, ^(He)center- 

9. Magnitude of the change in log(M/MQ) due to winds or accretion. 

10. Mass accreted so that compressional heating is correct (see §5.3). 

11. Changes in the logarithm of the total angular momentum. 

12. Distance moved in the HR diagram. 

13. Maximum allowed timestep under any circumstance. 

14. Any user specifed timestep limit, accomplished by setting raax_years_for .timestep, in the optional routine 
extras_checkjnodel. 

For convergence studies with respect to the timestep it is vital to change the control parameters that are actually setting the 
timestep. Often, this is just varcontrol_target, but in many situations the timestep will be set by one of the special timestep 
control parameters. 

Mesh Controls 

Control of the spatial mesh is a key ingredient of a stellar evolution instrument, and requires careful trade-offs. The mesh must 
respond to gradients in the structure, chemical composition, and energy generation, in order to give an accurate result, but it 
should not be overly dense since that will unnecessarily increase the cost of the calculation. 

Since MESA star allows for simulations with a fixed inner core mass, M^, the total mass M is M(. + where is the 
modeled mass. For cell k, MESA star stores the relative cell mass dqu - drnj^/M^ where dm^ is the mass contained in cell 
k. The relative mass interior to the outer cell face is - 1 - 2;=* ' d^/^ ^nd the total mass interior to the outer cell face is 
mil - qk * Mm + Mc. In all cases, m\ - M and qi = 1. We explicitly keep dq^ in addition to q and define q in terms of d^^ to 
avoid the need for evaluating q^ - qk+\ since that can involve the subtraction of almost equal numbers leading to an undesirable 
loss of precision (Lesaffre et al. 2006). For example, in the outer envelope of a star where the qu approach 1, the Aq^ can be 10"'^ 
or smaller. By storing dq^ we have 16 digit precision, whereas in this case, q^ - qk+i would only give us 4 digits at best for the 
relative cell mass. 

MESA star checks the structure and composition profiles of the model at the beginning of each timestep and, if necessary, 
adjusts the mesh. A single cell can be split into two or more cells and two or more adjacent cells can be merged. In practice, only 
a small fraction of the cells are changed during a remesh. This minimizes numerical diffusion, aids convergence, and keeps the 
cost of remeshing relatively small. Remeshing is divided into a planning stage and an adjustment stage. 

The planning stage determines which cells to split or merge based on the magnitude of allowed cell-to-cell changes in a variety 
of mesh functions. Built-in mesh functions include gradients of the mass, radius, pressure, temperature, adiabatic gradient, 
angular velocity and mass fractions above some threshold. Users can add others by defining their own other jnesh_functions 
routine. 

Other controls are provided to increase the sensitivity in regions selected by the user. Examples include increasing the spatial 
resolution in regions with changes in user-specified abundances with respect to pressure, changes in the energy generation rate 
with respect to pressure for different types of burning (e.g., the pp chains, CNO cycles, triple-or, and others), for regions near 
burning or non-burning convective boundaries, and others. 

After the mesh functions are evaluated, the relative magnitude of the changes between adjacent cells are determined. The 
magnitude of change is multiplied by mesh_delta_coe£f to obtain a weighted mesh function. Cells where the weighted changes 
are "too large" are marked for splitting, and cells where the changes are "too small" are marked for merging. For example, if the 
weighted changes in all mesh functions from cells k to k + n are less than 1, the series of cells from k to k + n are marked for 
merging. If any weighted mesh function changes from cell ktok+l by an amount greater than 1, the larger of cell k and cell k+ 1 is 
marked for splitting. Finally, if adjacent cells have too large of a relative size difference (as defined by meshjiiax_allowed_ratio 
which defaults to 2.5), the larger cell is marked for splitting and the check for excessive ratios is repeated. This can lead to a 
cascade of splitting in order to ensure that cells sizes do not have excessive jumps. 

The adjustment stage executes the remesh plan by performing the merge and split operations to calculate new values for 
basic variables. Special care is taken to use physical knowledge whenever possible when setting new values. For example, 
conservation of mass is accounted for when determining new densities, and species conservation is used when setting new mass 
fractions. Energy conservation is used when setting the temperature (see Paper I), and conservation of angular momentum plays 
a role in determining the angular velocity. Cells to be split are constructed by first performing a monotonicity-preserving cubic 
interpolation (Steffen 1990) in mass to obtain the luminosities and Inr values at the new cell boundaries. The new densities 
are then calculated from the new cell masses and volumes. Next, new composition mass fraction vectors are calculated. For 
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cells being merged, the mass averaged abundances are used. For cells being split, neighboring cells are used to form a linear 
approximation of mass fraction for each species as a function of mass coordinate within the cell. The slopes are adjusted so that 
the mass fractions sum to one everywhere, and the functions are integrated over the new cell mass to determine the abundances. 

Mass Adjustment 

Mass adjustment for mass loss or accretion is performed at each timestep when do_evolve_step calls the routine 
do.adjustjnass (see §B.l). MESA star offers several ways to set the rate of mass change M. A constant mass accretion 
rate (positive M) or mass loss rate (negative M) can be specified in the input files, a wind can produce a mass loss, the user can 
set M in an other_wind routine or in an other _checkjnodel routine. When do_adjust_mass is called, the timestep dt and the 
rate of mass change M are known, and thus the change in mass, 6M = M St. 

When there is a change in mass, instead of adding or removing cells, the total mass is changed by modifying the modeled 
mass Mn,, and cell mass sizes are changed by revisions to d^^ which in turn changes cell mass locations q/^ (see §B.4). The mass 
structure is divided into an inner region where the nik and dnit are unchanged but the qt and dqt change, an outer region where 
the qii and d^^ are unchanged but the and dm^ change, and an intermediate blending region where all of these change. The 
selection of the region boundaries is discussed in detail in Paper I. The implementation of e^av in the newly accreted matter is 
described in §5.3. 

Once the three regions have been defined, the dqi^ are updated. In the inner region they are rescaled by M/{M + 6M). Thus, dm^, 
111),, and X)^ have the same value before and after a change in mass. This eliminates the possibility of unwanted numerical diffusion 
causing unphysical mixing in the center region. In the outer region, cells retain the same value of dqi^ to improve convergence 
in the high entropy parts of the star (Sugimoto et al. 1981). In the intermediate region, the dqi^ are uniformly scaled to make 
Y,dqk = 1. 

The chemical mass fractions of cells in the intermediate and outer regions are then updated by summing the abundances 
between the new cell mass boundaries. This step is not necessary for the inner region since those cells have not changed mass 
location. In the case of mass accretion, the composition of the outermost cells whose enclosed mass totals 6M is set to match the 
specified accretion abundances. The single cell that is part old material and part newly accreted material is given an appropriately 
mixed composition. 

Finally, to create a somewhat better starting model for the Newton iterations (see §B.2), the InT and Inp and lni|as values 
are revised by monotonic cubic interpolating to the cell center by mass from the values prior to mass adjustment. The In r and 
material speed v are also set by monotonic cubic interpolation to the value at the new outer mass boundary. The angular velocity 
is set by integrating the angular momentum between the new cell mass boundaries and using the new In r values, conserving the 
total angular momentum to the floating point limit of the arithmetic. 

Evolving the Angular Velocity 

Initialization of rotation in MESA star begins from a non-rotating model. The angular velocity u is added to the set of model 
variables and initialized to a constant value throughout the model (i.e., solid body rotation). The initial value of <jj can be specified 
as a surface rotational velocity (in km/s) or as a fraction of the critical rotation rate (see §6). During the subsequent evolution, 
u) is changed at each timestep by remeshing, mass adjustment, radius adjustment (as part of the structure evolution), optional 
extra angular momentum removal in the outer layers, and the transport of angular momentum optionally with user-defined source 
terms for external torques. 

The angular velocity u is defined at cell boundaries. Thus omega (k) is at the outer boundary of cell k, which is the same 
location as the radius, r (k) , the specific moment of inertia, i _rot (k) , and the specific angular momentum, j _rot (k) . The mass 
associated with omega (k) spans the range from the center by mass of cell k outward to the center by mass of cell k-1 and is 
referred to as din_bar(k) to distinguish it from the cell mass dm(k). 

The remeshing operation splits and merges cells but does not change the physical stellar structure (see §B.4). For regions 
where there has been a change in the mesh, the values of w are adjusted to give the same angular momentum as before. More 
specifically, the angular momentum from the original model is summed over the mass range encompassed by the new dm_bar (k) , 
and omega (k) is adjusted to give the same total for the new model. 

During the mass adjustment operation, when mass is added or removed from the model, cells in the outer layers are moved to 
new mass locations (see §B.5). As part of this process, the angular velocity values are updated to conserve angular momentum 
using the same scheme as for remeshing: sum the angular momentum in the original model and set omega (k) in the new model 
to conserve it. Newly added material from accretion is given the current surface angular velocity. In the case of mass loss, this 
operation removes the amount of angular momentum contained in the lost mass at the start of the timestep; it does not deal with 
possible transport of angular momentum into the lost mass during the timestep. That is dealt with by an optional, user-specified 
removal prior to the angular momentum transport. 

MESA star performs the transport of angular momentum as a separate operation from the evolution of structure and composi- 
tion. This is done in order to obtain high accuracy in the angular momentum transport by using substeps and quad-precision linear 
algebra. It does not introduce additional operator splitting errors since u is not used in the structure and abundance equations. So 
we solve for the new structure and composition after any mass change and before the transport of angular momentum. Calculation 
of the new stellar structure changes the radii but does not change the mass partitioning of the model (see §B.2). Given the new 
radius r(k), we calculate the new i_rot(k). Then using the unchanged j_rot(k), omega(k) is set to j_rot(k)/i_rot(k) to 
conserve specific angular momentum. Since dm_bar (k) has not changed, this also conserves total angular momentum. 

Next, MESA star applies an optional, user-specified amount of angular momentum loss in the outer surface layers. This is to 
account for possible transport of angular momentum during the timestep from these outer layers into the mass removed by the 
mass adjustment operation. 
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The final operation is the transport of angular momentum within the star, which is treated with a diffusion approximation (Endal 
& Sofia 1978; Pinsonneault et al. 1989; Heger et al. 2000) 



dtj,,, i\dm 



(4;rrV) iv 



du) 
dm 



2ui3A /IdhW, ^^^^ 



r \dt /jj, \ 2 d In r 



where ; is the specific moment of inertia of a shell at mass coordinate m, and v is the turbulent viscosity determined as the sum 
of the diffusion coefficients for convection, double diffusion, overshooting and rotationally-induced instabilities (see §6). The 
diffusive transport is carefully implemented to accurately conserve angular momentum. The angular momentum associated with 
location k is dm_bar(k) "i_rot(k) "omega(k). The change in angular momentum for k is determined by the flux in angular 
momentum from A; - 1 to ^ and from k + \ io k. The flux from k - \ Xo k k set by v{k - 1) and the difference between 
omegaCk) and omega(k-l). The flux from + 1 to A: is found similarly using v{k) and the difference between omega(k) and 
omega (k+1). Source terms for location k are applied by user-supplied values for extra_jdot (k) or extra_oinegadot(k). The 
finite difference equation for the effects of the transport and source terms is solved over the stellar timestep with an implicit time 
integration that uses multiple smaller timesteps. The sizes of these substeps are determined by the timescale set by the diffusion 
coefficients and the differences in w. It is not unusual to use 10 or more substeps to evolve omega (k) over the stellar timestep. 
Each implicit substep is solved using a quad-precision tridiagonal matrix routine. The conservation of total angular momentum 
is monitored and the stellar timestep is rejected if there is any deviation from conservation by more than a user-specified factor. 
In practice, we find the total angular momentum is conserved over the stellar timestep to within a few digits of the floating point 
limit of the arithmetic. 



Nuclear Reactions 

A reaction network is defined by a set of isotopes and a set of reactions; these sets are specified in a reaction network 
definition file. MESA comes with many predefined reaction networks in data/net_data/nets and can also incorporate 
user-defined networks. To use a custom network, a user creates a reaction network definition file containing the command 
add_isos_and_reactions(isos_list), which will automatically add all reactions linking the isotopes in isos_list. The 
sequence of isotopes in isos_list may be specified by the name of the isotope: for example, add_isos_and_reactions(he4) 
adds "^He. Alternatively, one can specify the name of element followed by the desired minimum and maximum nucleon number. 
For example, the command add_isos_and_reactions(o 16 18) adds ""O, '^O, and "^O. Note that because many of the pre- 
defined networks may use effective rates — that is, using one reaction rate to represent a reaction sequence or group of reaction 
sequences — it is not recommended that the user extend one of the pre-existing networks with this command. 

MESA creates and stores reaction rate tables for each reaction whose entries are derived from evaluating standard analytic fitting 
formulas (see §A.4), but these reaction rates may be replaced with user-specified values. To change a rate for a given reaction, 

1. create a file with two columns: the temperature in units of 10** K and the rate Np^{crv) in units of cm^ g"' s"' ; 

2. list the file name in a local file rate_list.txt along with its "handle" for the reaction rate in question (see discussion 
below); and 

3. set the parameter rate_tables_dir in the namelist star_job to the name of the directory in which rate_list.txt is 
located; by default this is data/rates/rate_tables. 

The handle for a reaction is derived from the input and output channel isotopes according to a few rules. Capture reactions, such 
as x{p, y)y, have handles of the form r_x_pg_y and exchange reactions, such as x(a, p)y, have handles of the form r_x_ap_y. Other 
arbitrary reactions may be added by listing them in a form r_inputs_to_outputs where inputs and outputs are isotopes sepa- 
rated by If the same isotope appears two or more times, the isotope name may be repeated. For example, the triple-ff reaction 
is specified as r Jie4Jie4_he4_to_cl2. Isotopes are ordered by increasing Z and A^, e.g., r_h3_be7_tojieut Jil_he4 Jie4. To 
see a list of reactions used, the parameter show_net_reactions_info in namelist star, job should be set to . true . . 

Multicore Performance 

MESA implements shared memory multiprocessing via OpenMP". Paper I explored the runtime scaling of MESA star which 
at that time used a banded matrix linear algebra solver that did not benefit from multiple cores. A large part of the performance 
improvement in MESA star since Paper I comes from converting to a parallel block tridiagonal linear algebra solver derived from 
BCYCLIC (Hirshman et al. 2010). This improved solver is particularly important since linear algebra is typically the largest part 
of the runtime in MESA star. In addition, the new algorithm has the desirable property of producing numerically identical results 
independent of the number of cores, an attribute that is not generally true of parallel matrix solvers. 

Our test case is a 1.5 M© model with Z = 0.02 that is evolved from the ZAMS until the central H mass fraction falls to 0.35. 
This model includes 25 isotopes and 4 structure variables per cell with a variable number of zones typically exceeding 1700. The 
test takes ~55 time steps to cover ~ 1 .4 Gyr and uses the default amount of I/O. 

Figure 48 shows the scaling behavior of some key components of MESA star under GFORTRAN 4.7.2. on a 12 core 2010 
Apple MacPro. The dotted line shows the ideal scaling relation where doubling the number of cores cuts the run time in half. 
The linear algebra, labeled "mtx," dominates the total run time as the number of cores increases. For example, in the case of 12 
cores it accounts for about half the total and is 2.5 times larger than "net", the evaluation of the nuclear reaction network. The 

" http : //www . opermip . org 
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net evaluations closely approach the ideal scaling behavior because they can be done in parallel, each cell independent of the 
others, with one core working on one cell at a time. The equation of state component, labeled "eos," also closely approaches the 
ideal scaling law while consuming less than a third of the run time for the net. The component labeled "eqns," which includes 
the evaluation of the structure equations and the creation of the block tridiagonal matrix, also is close to the ideal scaling law 
and costs about the same as the eos. The "other" component is everything else. It is dominated by processes that currently are 
not efficient to parallelize because of the relatively large overhead for OpenMP operations. Consequently it remains at roughly 
a constant run time independent of the number of cores. When a significantly larger number of cores per processor becomes 
available, the larger operations in this category will have to be reworked or they will dominate the total run time. 
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Figure 48. Scaling beliavior of various components of MESA star using 1, 2, 3, 6, and 12 cores. The top curve shows the total run time, and the lower solid 
curves show the run times for the components of the total. The dotted line shows the ideal scaling relation. 

The run time also depends on the hardware, the quality of the compiled code, and the efficiency of the OpenMP implementation. 
For example, when the test case was rerun using IFORT 12.0.4 on the same 12-core machine, we found that GFORTRAN was 
6% faster that IFORT on one core but took 19% longer than IFORT on 12 cores. As another example, we ran the test case under 
GFORTRAN 4.7.2 on a 40-core server. While we obtained a speedup of 5.6 in going from 1 to 12 cores for the 12 core machine 
(see Figure 48), we find a speedup of only 4.8 on the 40-core server in going from 1 to 12 cores. Moreover, the speedup per core 
dropped steeply beyond 8-12 cores on the 40-core server, confirming the expectation that much work will be required to make 
full use of machines with many cores. 

Operating System and Compiler Considerations 

We next consider the range of values that are obtained when running on different operating systems and with different compilers 
and options. Such a comparison is important because computer hardware and software are heterogeneous and, under some 
circumstances, may give different answers. 

For this comparison we use two different systems, both with Intel Xeon processors, and either Linux (Cent OS 5.7) or Mac OS 
X (10.7.5). Both systems have Fortran compilers from GNU (version 4.7.2) and Intel (version 12). In addition to these operating 
system and compiler configurations, we consider effects of optimization and OpenMP multiprocessing. The default for compiling 
MESA is to enable OpenMP and use optimization at the -02 level. Here we compare results across the different operating systems, 
compilers, and without OpenMP or without optimization (i.e., -00). These comparisons were carried out with MESA revision 
4692. The test run is the solar calibration test_suite case with an absolute tolerance of 10"^. 

The following is a summary of our findings. On both operating systems MESA compiled by the GNU compiler with optimization 
and OpenMP support produces runs that are identical, to machine precision, with runs compiled by the same compiler but with 
either OpenMP or optimization disabled. The same is not true of the Intel compiler, for which we find that the calibrated values 
vary slightly, but within the specified tolerances, for runs with different numbers of threads. MESA running on the Linux system 
compiled with either compiler and the default options gives calibrated values that are compatible within the specified tolerances, 
but not identical to, those obtained using the same compiler and compiler options on the Mac OS X system. 

We found that the GNU and Intel compilers with optimization and OpenMP support produce somewhat different results. The 
difference is present on both Mac OS X and Linux operating systems. For our test configuration, MESA compiled with the Intel 
compiler produced values for the solar-calibrated model that are smaller than those obtained from MESA compiled with the GNU 
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Table 6 

Principal components of the MESA Software Development Kit 



Name 


Purpose 


Version 


License" 


GFORTRAN 


Compiler 


4.7.2 


Open source (GPL ver. 2) 


BLAS 


Matrix algebra 


2011-04-19 


Open source (other) 


LAPACK 


Matrix algebra 


3.4.2 


Open source (other) 


HDF5 


File storage 


1.8.9 


Open source (other) 


NDIFF 


Numerical comparison 


2.00 


Open source (GPL ver. 2) 


PGPLOT 


Plotting 


5.2.2 


Open source (non-commercial) 


SE 


File storage 


1.2.1 


Open source (other) 



""GPL" denotes the GNU General Public License (with the version in parentheses); "non-commercial" denotes an open-source license with restrictions on 
commercial distribution; and "other" denotes to a variety of open-source licenses which permit largely unrestricted distribution. 

compiler by ~ 0.01 in both initial Y and aMLT- Efforts will be made in the future to uncover the causes of these differences and 
eliminate them if possible. 

Visualization 

MESA star provides alphanumeric output at user-specified regular intervals. In addition, the routines in module 
star/public/pgstar . f provide an option for concurrent graphical output with the PGPLOT'^ library to create on-screen 
plots that can be saved for post-processing into animations of an evolutionary sequence. A variety of options are provided and 
are all configurable through the PGstar inlist. For example, a PGstar XI 1 window can simultaneously hold an H-R diagram, a 
Tc-pc diagram, and interior profiles of physical variables, such as nuclear energy generation and composition. The PGstar inlist 
is read at each timestep, so the display options can be changed without have to stop MESA star. 

Since Paper I, a number of MESA star users have developed and released toolkits'^ to visualize the alphanumeric output 
with common graphical packages including: Mathematica scripts (contributed by Richard O'Shaughnessy) and the intuitive 
and efficient graphical user interface MESAFace (Giannotti et al. 2012); MatLab utilities (contributed by Dave Spiegel and 
Gongjie Li); IDL functions (contributed by Rich Townsend); Python scripts (contributions from Falk Herwig and the NuGrid 
collaboration, David Kaplan, Alfred Gautschy, William Wolf); and Tioga scripts (contributed by Christopher Mankovich and Bill 
Pax ton). 

THE MESA SOFTWARE DEVELOPMENT KIT 

MESA is provided as source code, allowing users access to all of the implementation details. Installation necessarily involves 
building the code from source, which is a non-trivial task. A successful build requires cooperation between the operating system, 
compiler, libraries, and utilities. 

To address this issue we have created the MESA Software Development Kit (SDK), which packages everything necessary to 
establish a unified and maintained build environment. The principal components of the SDK are summarized in Table 6; all of these 
are distributed under an open-source license (detailed in the table), permitting their redistribution without financial or copyright 
encumbrances. Perhaps the most important component is the GFORTRAN compiler, part of the GNU Compiler Collection. 
GFORTRAN implements almost all of the Fortran 2003 (F2003) standard, and benefits from a high level of community support. 

The SDK is available for Intel x86 and x86-64 CPU architectures running the Linux and Mac OS X operating systems (these 
platforms comprise most of the MESA user base). Installation of the kit is straightforward, requiring a tar archive to be unpacked 
(Linux) or an application folder to be copied (OS X), followed by the initialization of a few environment variables. By default, 
MESA is configured to compile "out-of-the-box" with the SDK. MESA can also be compiled without the SDK, using any alternate 
compiler which supports the F2003 standard. In this respect, GFORTRAN should not be viewed as the MESA compiler (nor the 
full SDK as the MESA build environment). MESA will adhere to Fortran standards rather than rely on vendor-specific extensions. 

Uptake of the SDK has been very rapid: at the time of writing, we estimate over 90% of the MESA community (over 500 
users) are using the SDK. This growth has been matched by a significant decline in the number of installation support requests, 
and a corresponding reduction in the time taken to resolve these requests. With these maintenance overheads curbed, the MESA 
developers are able to devote more of their time to refining and extending the code. 
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